sn-ENS neuron integration

source("/Shared_win/projects/RNA_normal/analysis.10x.r")

load processed obj

GEX.seur <- readRDS("./sn10x_SZJ.sct_anno.rds")
GEX.seur
## An object of class Seurat 
## 47749 features across 28452 samples within 3 assays 
## Active assay: SCT (21008 features, 0 variable features)
##  2 other assays present: RNA, integrated
##  3 dimensional reductions calculated: pca, tsne, umap
DefaultAssay(GEX.seur) <- 'integrated'
DefaultAssay(GEX.seur)
## [1] "integrated"

filtering and re-clustering

GEX.seur <- subset(GEX.seur, subset= seurat_clusters %in% c(0:29))
GEX.seur <- subset(GEX.seur, subset= percent.mt <1.5)
GEX.seur <- subset(GEX.seur, subset= (cnt %in% c("Stst.CTL","Stst.CKO") & nFeature_RNA > 500 & nFeature_RNA < 1800 & nCount_RNA < 4000) |
                                     (cnt %in% c("CR7d.CTL","CR7d.CKO") & nFeature_RNA > 800 & nFeature_RNA < 2100 & nCount_RNA < 5000))
GEX.seur
## An object of class Seurat 
## 47749 features across 24358 samples within 3 assays 
## Active assay: integrated (2386 features, 2386 variable features)
##  2 other assays present: RNA, SCT
##  3 dimensional reductions calculated: pca, tsne, umap
length(GEX.seur@assays$integrated@var.features)
## [1] 2386
GEX.seur <- ScaleData(object = GEX.seur, verbose = TRUE, 
                         vars.to.regress = c("percent.mt","percent.rb","nCount_RNA"))
## Regressing out percent.mt, percent.rb, nCount_RNA
## Centering and scaling data matrix
GEX.seur <- RunPCA(GEX.seur, do.print = TRUE,
                      features = GEX.seur@assays$integrated@var.features,
                      seed.use = 133,
                      npcs = 100,
                      #ndims.print = 1,
                      verbose = T)
## PC_ 1 
## Positive:  Nrg3, Grid2, Lrrtm4, Pcdh15, Trpm3, Nrg3os, Ndst4, Cadm1, Kcnq3, Nkain2 
##     Galntl6, Cacnb2, Slit2, Ctnna3, Frmpd4, Epha5, Magi1, Syn3, Kcnma1, Fstl5 
##     Pde4b, Kcnc2, Plce1, Grip1, Tox, Kcnd2, Ryr2, Sntb1, Grik2, Lrp1b 
## Negative:  Ntrk3, Robo2, Tmeff2, Nrxn3, Ano2, Cdh8, Cpne4, Mgat4c, Myl1, Clstn2 
##     Plxna4, Pdzrn4, Zfp804a, Adgrg6, Spock3, Dgkg, Gpr149, Cdh6, Pcdh10, Cux2 
##     Astn2, Cntn5, Cysltr2, Kcnb2, Iqgap2, Arhgap6, Grin3a, Cacna2d3, Ccbe1, Itgb6 
## PC_ 2 
## Positive:  Rbfox1, Bnc2, Ptprt, Gpc6, Tafa1, Tshz2, Grik1, Mdga2, St6galnac3, Tox 
##     Adgrb3, Frmd4b, Brinp2, Fbxw15, Pcdh7, Cdc14a, Plcxd3, Agbl4, Fbxw24, Pbx1 
##     Unc5c, Chat, Pld5, Oprk1, Negr1, Dock2, Pde4b, Gfra2, Adamtsl1, Caln1 
## Negative:  Nos1, Auts2, Etv1, Cadps2, Gfra1, Alcam, Fam155a, Kcnab1, Egfem1, Kcnq5 
##     Asic2, Dgkb, Dach1, Plekha5, Schip1, Rgs6, Kcnt2, Epha5, Cmah, Ank2 
##     Nav3, Stxbp6, Cntnap5a, Creb5, Hs6st3, Tmem108, Lrrc4c, Ncam2, Frmpd4, Ablim2 
## PC_ 3 
## Positive:  Kcnip4, Cdh18, Csmd3, Kctd8, Cadm2, Klhl1, Gabrg3, Cntn3, Pbx3, Htr4 
##     Pde4d, Dlc1, Prkg1, Dmd, Dock10, Khdrbs2, Car10, Edil3, Meis1, Gpc6 
##     Skap1, Tenm2, Serpini1, Tac1, Plcl1, Grm7, C79798, Gda, March1, Nrp2 
## Negative:  Sgcd, Ptprg, Adgrg6, Nfia, Filip1, Cysltr2, Ccbe1, Slc4a4, Nfib, Gpr149 
##     Fgf13, Nmu, Sema6d, Ano2, Dgkg, Grin3a, Cpne4, Dapk2, Itgb6, Malat1 
##     Zfp804a, Ngfr, Cbln2, Efr3a, Nos1, Tshz2, Lhfpl2, Airn, Cntnap5a, Zfp521 
## PC_ 4 
## Positive:  Kcnt2, Lingo2, Fgf13, Prkg1, Dock10, Ndst4, Epha5, Lrrc4c, Ctnna3, Cntn5 
##     Lrrtm4, Dmd, Gda, Kcnip4, Tac1, Chl1, Ank2, Sorcs1, Rgs6, Nxph1 
##     Ptprz1, Hs3st2, Thsd4, Necab1, Grem2, Rora, Kctd8, Man1a, Galntl6, Sorcs3 
## Negative:  Trhde, Chsy3, Ebf1, Trps1, March1, Gal, Cntn4, Col18a1, Nrg1, Enox1 
##     Trpm3, Zmat4, Ntng1, Cpa6, Sdk1, Csmd1, Dcc, Shisa6, Ccser1, Tenm4 
##     Npas3, Nkain3, Kcnd2, Plcxd3, Tenm1, Sctr, Nwd2, Kcnh7, Prune2, Sez6l 
## PC_ 5 
## Positive:  Ptprd, Nrg1, Trhde, Lsamp, Egfem1, Cntn4, Cntn3, Adgrl2, Rmst, Cntn5 
##     Csmd3, Sgcz, Kcnd2, Luzp2, Gal, Nav2, Astn2, Ebf1, Asic2, Trps1 
##     Cpa6, Lrp1b, Kcnip4, Hs6st3, Sorcs1, Zbtb7c, Moxd1, Scn11a, Zmat4, Csmd1 
## Negative:  Dgkb, Klhl1, Vwc2l, Pbx3, Il1rapl1, Rasgef1b, Cdh9, Alk, Zfhx3, Mgat4c 
##     Sema5a, Dpp10, Vcan, Alcam, Galnt18, Galnt13, Auts2, Fam155a, Zbbx, C79798 
##     Olfr78, Thsd7b, Pcdh7, Scgn, Nek1, Serpini1, Sntg1, Lncbate1, P3h2, Mir100hg
DimHeatmap(GEX.seur, dims = 1:12, cells = 1500, balanced = TRUE,ncol = 4)

ElbowPlot(GEX.seur, ndims = 100)

ElbowPlot(GEX.seur, ndims = 50)

#
PCsct <- 1:18

GEX.seur <- FindNeighbors(GEX.seur, k.param = 20, dims = PCsct, compute.SNN = T, reduction = 'pca', verbose = T)
## Computing nearest neighbor graph
## Computing SNN
GEX.seur <- FindClusters(GEX.seur, dims.use = PCsct, algorithm = 1, save.SNN =T, resolution = 2, reduction = 'pca', verbose = T)
## Warning: The following arguments are not used: dims.use, save.SNN, reduction
## Suggested parameter: dims instead of dims.use
## Warning: The following arguments are not used: dims.use, save.SNN, reduction
## Suggested parameter: dims instead of dims.use
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 24358
## Number of edges: 917281
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8329
## Number of communities: 26
## Elapsed time: 4 seconds
GEX.seur <- RunTSNE(object = GEX.seur, assay = "integrated", seed.use = 233, dims = PCsct, complexity=100)
GEX.seur <- RunUMAP(object = GEX.seur, assay = "integrated", seed.use = 166, dims = PCsct, n.neighbors = 20, min.dist = 0.3)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 00:44:06 UMAP embedding parameters a = 0.9922 b = 1.112
## 00:44:06 Read 24358 rows and found 18 numeric columns
## 00:44:06 Using Annoy for neighbor search, n_neighbors = 20
## 00:44:06 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 00:44:08 Writing NN index file to temp file C:\Users\Shaorui\AppData\Local\Temp\RtmpMrRNx9\file848f084ce9
## 00:44:08 Searching Annoy index using 1 thread, search_k = 2000
## 00:44:14 Annoy recall = 100%
## 00:44:14 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 20
## 00:44:16 Initializing from normalized Laplacian + noise (using irlba)
## 00:44:18 Commencing optimization for 200 epochs, with 703756 positive edges
## 00:44:42 Optimization finished
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'tsne', group.by = "seurat_clusters") +
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters")

DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'tsne', group.by = "Anno1") +
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "Anno1")

DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'tsne', group.by = "Anno2") +
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "Anno2")

color.cnt <- scales::hue_pal()(4)[c(2,1,3,4)]
color.cnt
## [1] "#7CAE00" "#F8766D" "#00BFC4" "#C77CFF"
DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "cnt", split.by = "cnt", ncol = 4, cols = color.cnt)

DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "Anno1") +
DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "cnt",cols = color.cnt)

DimPlot(subset(GEX.seur, subset = cnt %in% c("Stst.CTL","Stst.CKO")), label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "cnt", cols = color.cnt[1:2]) +
DimPlot(subset(GEX.seur, subset = cnt %in% c("CR7d.CTL","CR7d.CKO")), label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "cnt", cols = color.cnt[3:4])

FeaturePlot(GEX.seur, 
            reduction = "umap", features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"))

FeaturePlot(subset(GEX.seur, subset = cnt %in% c("Stst.CTL","Stst.CKO")),
            reduction = "umap", features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"))

FeaturePlot(subset(GEX.seur, subset = cnt %in% c("CR7d.CTL","CR7d.CKO")), 
            reduction = "umap", features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"))

(DimPlot(subset(GEX.seur,subset=cnt %in% c("Stst.CTL","Stst.CKO")), label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "Anno1")+labs(title="Stst only")) +
(DimPlot(subset(GEX.seur,subset=cnt %in% c("CR7d.CTL","CR7d.CKO")), label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "Anno1")+labs(title="CR7d only"))

(DimPlot(subset(GEX.seur,subset=cnt %in% c("Stst.CTL","Stst.CKO")), label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters")+labs(title="Stst only")) +
(DimPlot(subset(GEX.seur,subset=cnt %in% c("CR7d.CTL","CR7d.CKO")), label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters")+labs(title="CR7d only"))

(DimPlot(subset(GEX.seur,subset=cnt %in% c("Stst.CTL","Stst.CKO")), label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "Anno1")+labs(title="Stst only")) +
(DimPlot(subset(GEX.seur,subset=cnt %in% c("Stst.CTL","Stst.CKO")), label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters")+labs(title="Stst only"))

(DimPlot(subset(GEX.seur,subset=cnt %in% c("CR7d.CTL","CR7d.CKO")), label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "Anno1")+labs(title="CR7d only")) +
(DimPlot(subset(GEX.seur,subset=cnt %in% c("CR7d.CTL","CR7d.CKO")), label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters")+labs(title="CR7d only"))

sort_clusters

GEX.seur$sort_clusters <- factor(as.character(GEX.seur$seurat_clusters),
                                 levels = c(3,0,7,8, 14, 13, 9, 1, 17,
                                            6,5,2,4, 12, 16, 23,
                                            10, 19,25, 22,
                                            18,11, 15,21, 24, 20))
VlnPlot(subset(GEX.seur, subset = cnt %in% c("Stst.CTL","Stst.CKO")), 
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "sort_clusters")

VlnPlot(subset(GEX.seur, subset = cnt %in% c("CR7d.CTL","CR7d.CKO")), 
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "sort_clusters")

previous markers

DefaultAssay(GEX.seur) <- "SCT"
DefaultAssay(GEX.seur)
## [1] "SCT"
markers.new.ss <- list(EMN=c("Chat","Bnc2","Tox","Ptprt",
                       "Gfra2","Oprk1","Adamtsl1", 
                       "Fbxw15","Fbxw24","Chrna7",
                       "Satb1","Itga6","Cntnap5b",
                       "Pgm5","Chgb",
                       "Nxph1",
                       "Gabrb1","Glp2r","Nebl",
                       "Lrrc7",
                       "Ryr3","Eda",
                       "Hgf","Lama2","Efnb2",
                       "Tac1",
                       "Kctd8",
                       "Ptn",
                       "Ntrk2","Penk","Itgb8",
                       "Fut9","Nfatc1","Egfr","Pdia5",
                       "Ahr","Mgll",
                       "Aff3",
                       "Chrm3"
                       ),
                 IMN=c("Nos1","Kcnab1",
                       "Gfra1","Etv1",
                       "Man1a","Airn",
                       "Adcy2",
                       "Col25a1",
                       "Cmah","Creb5","Vip","Pde1a",
                       "Ebf1","Gpc5","Mid1","Igfbp5",
                       "Ppara",
                       "Pcdh11x","Adcy8","Grp"
                       ),
                 IN=c("Npas3","Synpr","St18","Gal",
                      "Nova1",
                      "Cdh10","Kcnk13",
                      "Neurod6","Moxd1","Sctr",
                      "Piezo1","Vipr2","Adamts9","Sst","Kcnn2"
                      ),
                 IPAN=c("Calb2","Adcy1","Calcb","Nmu","Adgrg6","Pcdh10",
                        "Ngfr","Galr1","Il7","Aff2",
                        "Gpr149","Itgb6","Met","Itgbl1",
                        
                        "Cdh6","Cdh8",
                        "Clstn2","Ano2","Ntrk3",
                        "Cpne4","Vwc2l","Cdh9",
                        "Car10","Dcc",
                        "Scgn",
                        "Vcan","Cck","Piezo2","Kcnh7",
                        "Rerg","Bmpr1b","Skap1","Ntng1",
                        "Tafa2","Nxph2"),
                 CONTAM=c("Actl6b","Phox2b"), #,"Sox10","Plp1","L1cam","Gfap","Rxrg","Acta2","Actg2","Epcam","Pecam1","Ptprc"
                 pDEGs=c("Grid2","Ide","Btaf1","Slc15a2","Ccser1","Hdac9","Rspo2","Grem2"))

#
pm.All_new.s2 <- DotPlot(GEX.seur, features = as.vector(unlist(markers.new.ss)), group.by = "sort_clusters",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) #+ labs(title="All - newAnno1")
pm.All_new.s2

C25 ?

table(GEX.seur@meta.data[,c("sample","sort_clusters")])
##            sort_clusters
## sample        3   0   7   8  14  13   9   1  17   6   5   2   4  12  16  23  10
##   CR7d.CKO1 147  66 176 100  60  72  70 123  41 116  88 102 112  45  47  16 111
##   CR7d.CKO2  86  38 111  68  51  41  64  83  26  69  46  66  91  43  32  14  73
##   CR7d.CKO3 107  38 155 113  67  76  73 100  54 105  64 113 116  52  61  22  93
##   CR7d.CTL1 127  72 223 133  74  74  97 134  76 145  90 132 122  63  57  20  71
##   CR7d.CTL2 141  60 199 136  93  90  96 143  52 141  94 115 119  71  45  21  57
##   CR7d.CTL3 118  43 192 105  58  50  71 113  49 121  90 126 126  82  59  28  98
##   CR7d.CTL4 150  63 188 149  74  72 104 141  48 155  83 126 119  95  60  20  85
##   Stst.CKO1  95 244   4  45  69  68  93 135  53  58 132 146  92  88  48  44 129
##   Stst.CKO2 137 305   4  55  75 102  96 137  60  61 123 156 108  71  63  33 106
##   Stst.CKO3  99 224   2  64  66  79 110  99  54  62 113 110  76  62  61  22  93
##   Stst.CTL1 135 277   8  84  68  72 104 169  55 100 139 152  74 115  54  26  78
##   Stst.CTL2  93 174   6  71  46  55 100 111  59  68 106 112  89 102  42  24  69
##   Stst.CTL3 128 306  10  62  72  70  81 144  53  78 134 136  79  88  53  22  63
##            sort_clusters
## sample       19  25  22  18  11  15  21  24  20
##   CR7d.CKO1  42   2  22  47  75  41  45   8  49
##   CR7d.CKO2  30   6  23  34  70  40  29   9  20
##   CR7d.CKO3  42   9  37  31  75  52  37  11  40
##   CR7d.CTL1  48   4  19  41 101  62  41  15  42
##   CR7d.CTL2  29   5  31  23  72  49  51  21  34
##   CR7d.CTL3  60   5  28  40  67  45  60  10  36
##   CR7d.CTL4  46   4  42  48  97  67  48  22  48
##   Stst.CKO1  55   3  30  63  77  66  30  18  40
##   Stst.CKO2  51   5  17  65  98  80  25  20  51
##   Stst.CKO3  39   1  21  47  64  71  34  14  46
##   Stst.CTL1  51   1  15  44  65  66  35  15  42
##   Stst.CTL2  39   2  16  39  83  65  25  11  38
##   Stst.CTL3  39   3  16  67  90  70  35  17  57

I perfer to remove C25 which might be a rare new subtype, but here more like a contaminated mix in IN2

GEX.seur <- subset(GEX.seur, subset = seurat_clusters != c(25))
GEX.seur
## An object of class Seurat 
## 47749 features across 24308 samples within 3 assays 
## Active assay: SCT (21008 features, 0 variable features)
##  2 other assays present: RNA, integrated
##  3 dimensional reductions calculated: pca, tsne, umap
#
GEX.seur@meta.data[,grep("Anno|snn",colnames(GEX.seur@meta.data),value = T)] <- NULL
head(GEX.seur@meta.data)
##                        orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCAAGAATACAC-1_2 Stst.CTL_CKO       2282         1290  0.1314636  0.1314636
## AAACCCAAGCAATAGT-1_2 Stst.CTL_CKO       1535         1011  0.1302932  0.5863192
## AAACCCAAGGTGAGCT-1_2 Stst.CTL_CKO       1049          737  0.0000000  0.4766444
## AAACCCACAACGAGGT-1_2 Stst.CTL_CKO       1591         1046  0.6913891  0.4399749
## AAACCCACAAGAGTAT-1_2 Stst.CTL_CKO       2106         1250  0.2849003  0.2849003
## AAACCCACAATCAAGA-1_2 Stst.CTL_CKO       1451          953  0.1378360  0.1378360
##                      FB.info     S.Score    G2M.Score Phase      cnt  rep
## AAACCCAAGAATACAC-1_2   CTL.3 -0.00110357  0.012058898   G2M Stst.CTL rep3
## AAACCCAAGCAATAGT-1_2   CKO.3  0.02907907 -0.015558699     S Stst.CKO rep3
## AAACCCAAGGTGAGCT-1_2   CTL.1 -0.01100811 -0.008015087    G1 Stst.CTL rep1
## AAACCCACAACGAGGT-1_2   CTL.1 -0.02143685  0.005086498   G2M Stst.CTL rep1
## AAACCCACAAGAGTAT-1_2   CKO.3  0.02965845 -0.009057774     S Stst.CKO rep3
## AAACCCACAATCAAGA-1_2   CTL.2 -0.01216686  0.004143546   G2M Stst.CTL rep2
##                         sample tissue nCount_SCT nFeature_SCT seurat_clusters
## AAACCCAAGAATACAC-1_2 Stst.CTL3  Ileum       1860         1287               8
## AAACCCAAGCAATAGT-1_2 Stst.CKO3  Ileum       1525         1010               2
## AAACCCAAGGTGAGCT-1_2 Stst.CTL1  Ileum       1168          735               0
## AAACCCACAACGAGGT-1_2 Stst.CTL1  Ileum       1568         1045               3
## AAACCCACAAGAGTAT-1_2 Stst.CKO3  Ileum       1831         1248              16
## AAACCCACAATCAAGA-1_2 Stst.CTL2  Ileum       1451          953               3
##                      sort_clusters
## AAACCCAAGAATACAC-1_2             8
## AAACCCAAGCAATAGT-1_2             2
## AAACCCAAGGTGAGCT-1_2             0
## AAACCCACAACGAGGT-1_2             3
## AAACCCACAAGAGTAT-1_2            16
## AAACCCACAATCAAGA-1_2             3

intAnno

# intAnno1 
GEX.seur$intAnno1 <- as.character(GEX.seur$seurat_clusters)

GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(3,0,7,8,14)] <- "EMN1"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(13)] <- "EMN2"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(9)] <- "EMN3"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(1)] <- "EMN4"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(17)] <- "EMN5"

GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(6,5,2,4)] <- "IMN1"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(12)] <- "IMN2"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(16)] <- "IMN3"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(23)] <- "IMN4"

GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(10)] <- "IN1"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(19)] <- "IN2"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(22)] <- "IN3"

GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(18,11)] <- "IPAN1"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(15,21)] <- "IPAN2"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(24)] <- "IPAN3"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(20)] <- "IPAN4"

GEX.seur$intAnno1 <- factor(GEX.seur$intAnno1,
                            levels = c(paste0("EMN",1:5),
                                       paste0("IMN",1:4),
                                       paste0("IN",1:3),
                                       paste0("IPAN",1:4)))

# intAnno2
GEX.seur$intAnno2 <- as.character(GEX.seur$seurat_clusters)

GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(3,0,7,8,14)] <- "EMN1"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(13)] <- "EMN2"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(9)] <- "EMN3"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(1)] <- "EMN4"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(17)] <- "EMN5"

GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(6,5,2,4)] <- "IMN1"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(12)] <- "IMN2"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(16)] <- "IMN3"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(23)] <- "IMN4"

GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(10)] <- "IN1"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(19)] <- "IN2"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(22)] <- "IN3"

GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(18)] <- "IPAN1.1"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(11)] <- "IPAN1.2"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(15)] <- "IPAN2.1"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(21)] <- "IPAN2.2"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(24)] <- "IPAN3"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(20)] <- "IPAN4"

GEX.seur$intAnno2 <- factor(GEX.seur$intAnno2,
                            levels = c(paste0("EMN",1:5),
                                       paste0("IMN",1:4),
                                       paste0("IN",1:3),
                                       paste0("IPAN",c(1.1,1.2,2.1,2.2,3,4))))
#saveRDS(GEX.seur, "./sn10x_SZJ.sct_anno.s.rds")
#GEX.seur <- readRDS("./sn10x_SZJ.sct_anno.s.rds")
pp.umap1 <- DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'tsne', group.by = "seurat_clusters") +
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters")
pp.umap1

pp.umap2 <- DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "intAnno1") +
DimPlot(GEX.seur, label = T, label.size = 2.65, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "intAnno2")
pp.umap2 

pp.umap3 <- DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "cnt", split.by = "cnt", ncol = 4, cols = color.cnt)
pp.umap3

dotplot

markers.new.ss <- list(EMN=c("Chat","Bnc2","Tox","Ptprt",
                       "Gfra2","Oprk1","Adamtsl1", 
                       "Fbxw15","Fbxw24","Chrna7",
                       "Satb1","Itga6","Cntnap5b",
                       "Pgm5","Chgb",
                       "Nxph1",
                       "Gabrb1","Glp2r","Nebl",
                       "Lrrc7",
                       "Ryr3","Eda",
                       "Hgf","Lama2","Efnb2",
                       "Tac1",
                       "Kctd8",
                       "Ptn",
                       "Ntrk2","Penk","Itgb8",
                       "Fut9","Nfatc1","Egfr","Pdia5",
                       "Ahr","Mgll",
                       "Aff3",
                       "Chrm3"
                       ),
                 IMN=c("Nos1","Kcnab1",
                       "Gfra1","Etv1",
                       "Man1a","Airn",
                       "Adcy2",
                       "Col25a1",
                       "Cmah","Creb5","Vip","Pde1a",
                       "Ebf1","Gpc5","Mid1","Igfbp5",
                       "Ppara",
                       "Pcdh11x","Adcy8","Grp"
                       ),
                 IN=c("Npas3","Synpr","St18","Gal",
                      "Nova1",
                      "Cdh10","Kcnk13",
                      "Neurod6","Moxd1","Sctr",
                      "Piezo1","Vipr2","Adamts9","Sst","Kcnn2"
                      ),
                 IPAN=c("Calb2","Adcy1","Calcb","Nmu","Adgrg6","Pcdh10",
                        "Ngfr","Galr1","Il7","Aff2",
                        "Gpr149","Itgb6","Met","Itgbl1",
                        
                        "Cdh6","Cdh8",
                        "Clstn2","Ano2","Ntrk3",
                        "Cpne4","Vwc2l","Cdh9",
                        "Car10","Dcc",
                        "Scgn",
                        "Vcan","Cck","Piezo2","Kcnh7",
                        "Rerg","Bmpr1b","Skap1","Ntng1",
                        "Tafa2","Nxph2"),
                 CONTAM=c("Actl6b","Phox2b"),
                 pDEGs=c("Grid2","Ide","Btaf1","Slc15a2","Ccser1","Hdac9","Rspo2","Grem2"))
check.markers.ss <- as.vector(unlist(markers.new.ss))
check.markers.ss
##   [1] "Chat"     "Bnc2"     "Tox"      "Ptprt"    "Gfra2"    "Oprk1"   
##   [7] "Adamtsl1" "Fbxw15"   "Fbxw24"   "Chrna7"   "Satb1"    "Itga6"   
##  [13] "Cntnap5b" "Pgm5"     "Chgb"     "Nxph1"    "Gabrb1"   "Glp2r"   
##  [19] "Nebl"     "Lrrc7"    "Ryr3"     "Eda"      "Hgf"      "Lama2"   
##  [25] "Efnb2"    "Tac1"     "Kctd8"    "Ptn"      "Ntrk2"    "Penk"    
##  [31] "Itgb8"    "Fut9"     "Nfatc1"   "Egfr"     "Pdia5"    "Ahr"     
##  [37] "Mgll"     "Aff3"     "Chrm3"    "Nos1"     "Kcnab1"   "Gfra1"   
##  [43] "Etv1"     "Man1a"    "Airn"     "Adcy2"    "Col25a1"  "Cmah"    
##  [49] "Creb5"    "Vip"      "Pde1a"    "Ebf1"     "Gpc5"     "Mid1"    
##  [55] "Igfbp5"   "Ppara"    "Pcdh11x"  "Adcy8"    "Grp"      "Npas3"   
##  [61] "Synpr"    "St18"     "Gal"      "Nova1"    "Cdh10"    "Kcnk13"  
##  [67] "Neurod6"  "Moxd1"    "Sctr"     "Piezo1"   "Vipr2"    "Adamts9" 
##  [73] "Sst"      "Kcnn2"    "Calb2"    "Adcy1"    "Calcb"    "Nmu"     
##  [79] "Adgrg6"   "Pcdh10"   "Ngfr"     "Galr1"    "Il7"      "Aff2"    
##  [85] "Gpr149"   "Itgb6"    "Met"      "Itgbl1"   "Cdh6"     "Cdh8"    
##  [91] "Clstn2"   "Ano2"     "Ntrk3"    "Cpne4"    "Vwc2l"    "Cdh9"    
##  [97] "Car10"    "Dcc"      "Scgn"     "Vcan"     "Cck"      "Piezo2"  
## [103] "Kcnh7"    "Rerg"     "Bmpr1b"   "Skap1"    "Ntng1"    "Tafa2"   
## [109] "Nxph2"    "Actl6b"   "Phox2b"   "Grid2"    "Ide"      "Btaf1"   
## [115] "Slc15a2"  "Ccser1"   "Hdac9"    "Rspo2"    "Grem2"
length(check.markers.ss)
## [1] 119
sum(check.markers.ss %in% rownames(GEX.seur))
## [1] 119
pm.All_new.a <- DotPlot(GEX.seur, features = check.markers.ss, group.by = "sort_clusters",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All - sort_clusters")
pm.All_new.a

pm.All_new.b <- DotPlot(GEX.seur, features = check.markers.ss, group.by = "intAnno1",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All - intAnno1")
pm.All_new.b

pm.All_new.c <- DotPlot(GEX.seur, features = check.markers.ss, group.by = "intAnno2",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All - intAnno2")
pm.All_new.c

#
pm.All_new.c1 <- DotPlot(subset(GEX.seur, subset= cnt %in% c("Stst.CTL")), features = check.markers.ss, group.by = "intAnno2",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="Stst.CTL - intAnno2")
pm.All_new.c1

#
pm.All_new.c2 <- DotPlot(subset(GEX.seur, subset= cnt %in% c("Stst.CKO")), features = check.markers.ss, group.by = "intAnno2",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="Stst.CKO - intAnno2")
pm.All_new.c2

#
pm.All_new.c3 <- DotPlot(subset(GEX.seur, subset= cnt %in% c("CR7d.CTL")), features = check.markers.ss, group.by = "intAnno2",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="CR7d.CTL - intAnno2")
pm.All_new.c3

#
pm.All_new.c4 <- DotPlot(subset(GEX.seur, subset= cnt %in% c("CR7d.CKO")), features = check.markers.ss, group.by = "intAnno2",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="CR7d.CKO - intAnno2")
pm.All_new.c4

new

# define intAnno1/2 colors
color.A1 <- c("#678BB1","#8AB6CE","#3975C1","#669FDF","#4CC1BD",
              "#BF7E6B","#D46B35","#F19258","#FF8080",
              "#BDAE8D","#BD66C4","#C03778",
              "#97BA59","#DFAB16","#2BA956","#9FE727")
names(color.A1) <- levels(GEX.seur$intAnno1)

color.A2 <- c("#678BB1","#8AB6CE","#3975C1","#669FDF","#4CC1BD",
              "#BF7E6B","#D46B35","#F19258","#FF8080",
              "#BDAE8D","#BD66C4","#C03778",
              "#97BA59","#C4D116", "#DFAB16","#EDE25A", "#2BA956","#9FE727")
names(color.A2) <- levels(GEX.seur$intAnno2)

# define batch/condition colors
color.cnt <- scales::hue_pal()(4)[c(2,1,3,4)]
names(color.cnt) <- levels(GEX.seur$cnt)

color.test1 <- color.cnt[1:2]
color.test2 <- color.cnt[3:4]

## define feature colors
# Cell2020     
material.heat = function(n){
  colorRampPalette(
    c(
      "#283593", # indigo 800
      "#3F51B5", # indigo
      "#2196F3", # blue
      "#00BCD4", # cyan
      "#4CAF50", # green
      "#8BC34A", # light green
      "#CDDC39", # lime
      "#FFEB3B", # yellow
      "#FFC107", # amber
      "#FF9800", # orange
      "#FF5722"  # deep orange
    )
  )(n)
}

# Immunity2019, na gray
colors.Immunity <-c("#191970","#121285","#0C0C9A","#0707B0","#0101C5","#0014CF","#0033D3","#0053D8","#0072DD","#0092E1","#00B2E6",
                  "#00D1EB","#23E8CD","#7AF17B","#D2FA29","#FFEB00","#FFC300","#FF9B00","#FF8400","#FF7800","#FF6B00","#FF5F00","#FF5300",
                  "#FF4700","#F73B00","#EF2E00","#E62300","#DD1700","#D50B00","#CD0000")


# NatNeur2021, sc-neurons
color.ref <- c("#8AB6CE","#678BB1","#3975C1","#4CC1BD",
               "#C03778","#97BA59","#DFAB16","#BF7E6B",
               "#D46B35","#BDAE8D","#BD66C4","#2BA956",
               "#FF8080","#FF8080","#FF8080","#FF0000")
cowplot::plot_grid(
DimPlot(GEX.seur, reduction = "umap", group.by = "intAnno1", label = T, label.size = 3.25,repel = F, pt.size = 0.05,
        cols = color.A1),
  DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "cnt",
        cols = color.cnt) ,
rel_widths = c(4.8,5),
ncol = 2)

cowplot::plot_grid(
DimPlot(GEX.seur, reduction = "umap", group.by = "intAnno2", label = T, label.size = 2.65,repel = F, pt.size = 0.05,
        cols = color.A2),
  DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "cnt",
        cols = color.cnt) ,
rel_widths = c(4.85,5),
ncol = 2)

condition

Stst

test1.seur <- subset(GEX.seur, subset= cnt %in% c("Stst.CTL","Stst.CKO"))
test1.seur
## An object of class Seurat 
## 47749 features across 11409 samples within 3 assays 
## Active assay: SCT (21008 features, 0 variable features)
##  2 other assays present: RNA, integrated
##  3 dimensional reductions calculated: pca, tsne, umap
cowplot::plot_grid(
DimPlot(test1.seur, reduction = "umap", group.by = "intAnno1", label = T, label.size = 3.25,repel = F, pt.size = 0.15,
        cols = color.A1),
  DimPlot(test1.seur, label = F, pt.size = 0.15, repel = F, reduction = 'umap', group.by = "cnt",
        cols = color.test1) ,
rel_widths = c(4.8,5),
ncol = 2)

cowplot::plot_grid(
DimPlot(test1.seur, reduction = "umap", group.by = "intAnno2", label = T, label.size = 2.65,repel = F, pt.size = 0.15,
        cols = color.A2),
  DimPlot(test1.seur, label = F, pt.size = 0.15, repel = F, reduction = 'umap', group.by = "cnt",
        cols = color.test1) ,
rel_widths = c(4.85,5),
ncol = 2)

CR7d

test2.seur <- subset(GEX.seur, subset= cnt %in% c("CR7d.CTL","CR7d.CKO"))
test2.seur
## An object of class Seurat 
## 47749 features across 12899 samples within 3 assays 
## Active assay: SCT (21008 features, 0 variable features)
##  2 other assays present: RNA, integrated
##  3 dimensional reductions calculated: pca, tsne, umap
cowplot::plot_grid(
DimPlot(test2.seur, reduction = "umap", group.by = "intAnno1", label = T, label.size = 3.25,repel = F, pt.size = 0.15,
        cols = color.A1),
  DimPlot(test2.seur, label = F, pt.size = 0.15, repel = F, reduction = 'umap', group.by = "cnt",
        cols = color.test2) ,
rel_widths = c(4.8,5),
ncol = 2)

cowplot::plot_grid(
DimPlot(test2.seur, reduction = "umap", group.by = "intAnno2", label = T, label.size = 2.65,repel = F, pt.size = 0.15,
        cols = color.A2),
  DimPlot(test2.seur, label = F, pt.size = 0.15, repel = F, reduction = 'umap', group.by = "cnt",
        cols = color.test2) ,
rel_widths = c(4.85,5),
ncol = 2)

cell composition

Stst

cowplot::plot_grid(
pheatmap::pheatmap(table(cnt=test1.seur$sample,
      clusters=test1.seur$intAnno1)[c(4:6,1:3),],
                   main = "Cell Count",
                   gaps_row = c(3),
      gaps_col = c(5,9,12),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,

pheatmap::pheatmap(100*rowRatio(table(cnt=test1.seur$sample,
                                clusters=test1.seur$intAnno1)[c(4:6,1:3),] ),
                   main = "Cell Ratio",
                   gaps_row = c(3),
      gaps_col = c(5,9,12),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)

cowplot::plot_grid(
pheatmap::pheatmap(table(cnt=test1.seur$sample,
      clusters=test1.seur$intAnno2)[c(4:6,1:3),],
                   main = "Cell Count",
                   gaps_row = c(3),
      gaps_col = c(5,9,12),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,

pheatmap::pheatmap(100*rowRatio(table(cnt=test1.seur$sample,
                                clusters=test1.seur$intAnno2)[c(4:6,1:3),] ),
                   main = "Cell Ratio",
                   gaps_row = c(3),
      gaps_col = c(5,9,12),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)

barplot

stat1_intAnno1 <- test1.seur@meta.data[,c("cnt","rep","intAnno1")]

stat1_intAnno1.s <- stat1_intAnno1 %>%
  group_by(cnt,rep,intAnno1) %>%
  summarise(count=n()) %>%
  mutate(prop= count/sum(count)) %>%
  ungroup
## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.
pcom1.intAnno1 <- stat1_intAnno1.s %>%
  ggplot(aes(x = intAnno1, y = 100*prop, color=cnt)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.test1, name="") +
  labs(title="Cell Composition of Stst.CTL_CKO, intAnno1", y = "Proportion") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
  
  geom_point(aes(x=intAnno1, y = 100*prop, color= cnt),
             position = position_dodge(0.75),
             shape=16,alpha=0.75,size=2.15,
             stroke=0.15, show.legend=F)

pcom1.intAnno1

stat1_intAnno2 <- test1.seur@meta.data[,c("cnt","rep","intAnno2")]

stat1_intAnno2.s <- stat1_intAnno2 %>%
  group_by(cnt,rep,intAnno2) %>%
  summarise(count=n()) %>%
  mutate(prop= count/sum(count)) %>%
  ungroup
## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.
pcom1.intAnno2 <- stat1_intAnno2.s %>%
  ggplot(aes(x = intAnno2, y = 100*prop, color=cnt)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.test1, name="") +
  labs(title="Cell Composition of Stst.CTL_CKO, intAnno2", y = "Proportion") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
  
  geom_point(aes(x=intAnno2, y = 100*prop, color= cnt),
             position = position_dodge(0.75),
             shape=16,alpha=0.75,size=2.15,
             stroke=0.15, show.legend=F)

pcom1.intAnno2

sig.test

glm.nb - anova.LRT

Sort.N <- c("Stst.CTL","Stst.CKO")

glm.nb_anova.1.intAnno1 <- list()

for(x1 in 1:2){
  for(x2 in 1:2){
    if(x2>x1){
      
      stat1_intAnno1.s_N <- stat1_intAnno1.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
      
      total.N <- stat1_intAnno1.s_N %>%
        group_by(cnt, rep) %>%
        summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
      
      rownames(total.N) <- paste0(as.character(total.N$cnt),
                                  "_",
                                  as.character(total.N$rep))
      
      stat1_intAnno1.s_N$total <- total.N[paste0(stat1_intAnno1.s_N$cnt,
                                            "_",
                                            stat1_intAnno1.s_N$rep),"total"]
      
      glm.nb_anova.1.intAnno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
      
      for(AA in levels(stat1_intAnno1.s_N$intAnno1)){
        glm.nb_anova.1.intAnno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
          aaa <- stat1_intAnno1.s_N %>% filter(intAnno1==AA)
          aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
          aaa.anova <- anova(aaa.glmnb, test = "LRT")
          aaa.anova$'Pr(>Chi)'[2]
        }, error=function(e){
        
          tryCatch({
            aaa <- stat1_intAnno1.s_N %>% filter(intAnno1==AA)
            aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
            aaa.anova <- anova(aaa.glmnb, test = "LRT")
            aaa.anova$'Pr(>Chi)'[2]
          }, error=function(e){
            
            tryCatch({
              aaa <- stat1_intAnno1.s_N %>% filter(intAnno1==AA)
              aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
              aaa.anova <- anova(aaa.glmnb, test = "LRT")
              aaa.anova$'Pr(>Chi)'[2]
            }, error = function(e){
              NA
            })     
          })
        })
      }
      glm.nb_anova.1.intAnno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.1.intAnno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]]) 
    }
  }
}
glm.nb_anova.1.intAnno1_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.1.intAnno1)))
rownames(glm.nb_anova.1.intAnno1_df) <- names(glm.nb_anova.1.intAnno1)
colnames(glm.nb_anova.1.intAnno1_df) <- gsub("X","C",colnames(glm.nb_anova.1.intAnno1_df))
glm.nb_anova.1.intAnno1_df
##                        EMN1       EMN2      EMN3       EMN4      EMN5      IMN1
## Stst.CTLvsStst.CKO 0.437101 0.02226962 0.8008005 0.03400082 0.8765206 0.3068567
##                           IMN2      IMN3       IMN4          IN1       IN2
## Stst.CTLvsStst.CKO 0.002090672 0.2576225 0.07350664 2.335817e-05 0.4113903
##                           IN3     IPAN1     IPAN2     IPAN3    IPAN4
## Stst.CTLvsStst.CKO 0.06840589 0.7111369 0.8422819 0.4009551 0.888078
#options (scipen = 999)
round(glm.nb_anova.1.intAnno1_df,6)
##                        EMN1    EMN2     EMN3     EMN4     EMN5     IMN1
## Stst.CTLvsStst.CKO 0.437101 0.02227 0.800801 0.034001 0.876521 0.306857
##                        IMN2     IMN3     IMN4     IN1     IN2      IN3    IPAN1
## Stst.CTLvsStst.CKO 0.002091 0.257623 0.073507 2.3e-05 0.41139 0.068406 0.711137
##                       IPAN2    IPAN3    IPAN4
## Stst.CTLvsStst.CKO 0.842282 0.400955 0.888078
Sort.N <- c("Stst.CTL","Stst.CKO")

glm.nb_anova.1.intAnno2 <- list()

for(x1 in 1:2){
  for(x2 in 1:2){
    if(x2>x1){
      
      stat1_intAnno2.s_N <- stat1_intAnno2.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
      
      total.N <- stat1_intAnno2.s_N %>%
        group_by(cnt, rep) %>%
        summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
      
      rownames(total.N) <- paste0(as.character(total.N$cnt),
                                  "_",
                                  as.character(total.N$rep))
      
      stat1_intAnno2.s_N$total <- total.N[paste0(stat1_intAnno2.s_N$cnt,
                                            "_",
                                            stat1_intAnno2.s_N$rep),"total"]
      
      glm.nb_anova.1.intAnno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
      
      for(AA in levels(stat1_intAnno2.s_N$intAnno2)){
        glm.nb_anova.1.intAnno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
          aaa <- stat1_intAnno2.s_N %>% filter(intAnno2==AA)
          aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
          aaa.anova <- anova(aaa.glmnb, test = "LRT")
          aaa.anova$'Pr(>Chi)'[2]
        }, error=function(e){
        
          tryCatch({
            aaa <- stat1_intAnno2.s_N %>% filter(intAnno2==AA)
            aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
            aaa.anova <- anova(aaa.glmnb, test = "LRT")
            aaa.anova$'Pr(>Chi)'[2]
          }, error=function(e){
            
            tryCatch({
              aaa <- stat1_intAnno2.s_N %>% filter(intAnno2==AA)
              aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
              aaa.anova <- anova(aaa.glmnb, test = "LRT")
              aaa.anova$'Pr(>Chi)'[2]
            }, error = function(e){
              NA
            })     
          })
        })
      }
      glm.nb_anova.1.intAnno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.1.intAnno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]]) 
    }
  }
}
glm.nb_anova.1.intAnno2_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.1.intAnno2)))
rownames(glm.nb_anova.1.intAnno2_df) <- names(glm.nb_anova.1.intAnno2)
colnames(glm.nb_anova.1.intAnno2_df) <- gsub("X","C",colnames(glm.nb_anova.1.intAnno2_df))
glm.nb_anova.1.intAnno2_df
##                        EMN1       EMN2      EMN3       EMN4      EMN5      IMN1
## Stst.CTLvsStst.CKO 0.437101 0.02226962 0.8008005 0.03400082 0.8765206 0.3068567
##                           IMN2      IMN3       IMN4          IN1       IN2
## Stst.CTLvsStst.CKO 0.002090672 0.2576225 0.07350664 2.335817e-05 0.4113903
##                           IN3   IPAN1.1   IPAN1.2   IPAN2.1   IPAN2.2     IPAN3
## Stst.CTLvsStst.CKO 0.06840589 0.2661605 0.8564025 0.5426043 0.5770597 0.4009551
##                       IPAN4
## Stst.CTLvsStst.CKO 0.888078
#options (scipen = 999)
round(glm.nb_anova.1.intAnno2_df,5)
##                      EMN1    EMN2   EMN3  EMN4    EMN5    IMN1    IMN2    IMN3
## Stst.CTLvsStst.CKO 0.4371 0.02227 0.8008 0.034 0.87652 0.30686 0.00209 0.25762
##                       IMN4   IN1     IN2     IN3 IPAN1.1 IPAN1.2 IPAN2.1
## Stst.CTLvsStst.CKO 0.07351 2e-05 0.41139 0.06841 0.26616  0.8564  0.5426
##                    IPAN2.2   IPAN3   IPAN4
## Stst.CTLvsStst.CKO 0.57706 0.40096 0.88808

CR7d

heatmap

cowplot::plot_grid(
pheatmap::pheatmap(table(cnt=test2.seur$sample,
      clusters=test2.seur$intAnno1)[c(4:7,1:3),],
                   main = "Cell Count",
                   gaps_row = c(4),
      gaps_col = c(5,9,12),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,

pheatmap::pheatmap(100*rowRatio(table(cnt=test2.seur$sample,
                                clusters=test2.seur$intAnno1)[c(4:7,1:3),] ),
                   main = "Cell Ratio",
                   gaps_row = c(4),
      gaps_col = c(5,9,12),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)

cowplot::plot_grid(
pheatmap::pheatmap(table(cnt=test2.seur$sample,
      clusters=test2.seur$intAnno2)[c(4:7,1:3),],
                   main = "Cell Count",
                   gaps_row = c(4),
      gaps_col = c(5,9,12),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,

pheatmap::pheatmap(100*rowRatio(table(cnt=test2.seur$sample,
                                clusters=test2.seur$intAnno2)[c(4:7,1:3),] ),
                   main = "Cell Ratio",
                   gaps_row = c(4),
      gaps_col = c(5,9,12),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)

barplot

stat2_intAnno1 <- test2.seur@meta.data[,c("cnt","rep","intAnno1")]

stat2_intAnno1.s <- stat2_intAnno1 %>%
  group_by(cnt,rep,intAnno1) %>%
  summarise(count=n()) %>%
  mutate(prop= count/sum(count)) %>%
  ungroup
## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.
pcom2.intAnno1 <- stat2_intAnno1.s %>%
  ggplot(aes(x = intAnno1, y = 100*prop, color=cnt)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.test2, name="") +
  labs(title="Cell Composition of CR7d.CTL_CKO, intAnno1", y = "Proportion") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
  
  geom_point(aes(x=intAnno1, y = 100*prop, color= cnt),
             position = position_dodge(0.75),
             shape=16,alpha=0.75,size=2.15,
             stroke=0.15, show.legend=F)

pcom2.intAnno1

stat2_intAnno2 <- test2.seur@meta.data[,c("cnt","rep","intAnno2")]

stat2_intAnno2.s <- stat2_intAnno2 %>%
  group_by(cnt,rep,intAnno2) %>%
  summarise(count=n()) %>%
  mutate(prop= count/sum(count)) %>%
  ungroup
## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.
pcom2.intAnno2 <- stat2_intAnno2.s %>%
  ggplot(aes(x = intAnno2, y = 100*prop, color=cnt)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.test2, name="") +
  labs(title="Cell Composition of CR7d.CTL_CKO, intAnno2", y = "Proportion") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
  
  geom_point(aes(x=intAnno2, y = 100*prop, color= cnt),
             position = position_dodge(0.75),
             shape=16,alpha=0.75,size=2.15,
             stroke=0.15, show.legend=F)

pcom2.intAnno2

sig.test

glm.nb - anova.LRT

Sort.N <- c("CR7d.CTL","CR7d.CKO")

glm.nb_anova.2.intAnno1 <- list()

for(x1 in 1:2){
  for(x2 in 1:2){
    if(x2>x1){
      
      stat2_intAnno1.s_N <- stat2_intAnno1.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
      
      total.N <- stat2_intAnno1.s_N %>%
        group_by(cnt, rep) %>%
        summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
      
      rownames(total.N) <- paste0(as.character(total.N$cnt),
                                  "_",
                                  as.character(total.N$rep))
      
      stat2_intAnno1.s_N$total <- total.N[paste0(stat2_intAnno1.s_N$cnt,
                                            "_",
                                            stat2_intAnno1.s_N$rep),"total"]
      
      glm.nb_anova.2.intAnno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
      
      for(AA in levels(stat2_intAnno1.s_N$intAnno1)){
        glm.nb_anova.2.intAnno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
          aaa <- stat2_intAnno1.s_N %>% filter(intAnno1==AA)
          aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
          aaa.anova <- anova(aaa.glmnb, test = "LRT")
          aaa.anova$'Pr(>Chi)'[2]
        }, error=function(e){
        
          tryCatch({
            aaa <- stat2_intAnno1.s_N %>% filter(intAnno1==AA)
            aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
            aaa.anova <- anova(aaa.glmnb, test = "LRT")
            aaa.anova$'Pr(>Chi)'[2]
          }, error=function(e){
            
            tryCatch({
              aaa <- stat2_intAnno1.s_N %>% filter(intAnno1==AA)
              aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
              aaa.anova <- anova(aaa.glmnb, test = "LRT")
              aaa.anova$'Pr(>Chi)'[2]
            }, error = function(e){
              NA
            })     
          })
        })
      }
      glm.nb_anova.2.intAnno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.2.intAnno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]]) 
    }
  }
}
glm.nb_anova.2.intAnno1_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.2.intAnno1)))
rownames(glm.nb_anova.2.intAnno1_df) <- names(glm.nb_anova.2.intAnno1)
colnames(glm.nb_anova.2.intAnno1_df) <- gsub("X","C",colnames(glm.nb_anova.2.intAnno1_df))
glm.nb_anova.2.intAnno1_df
##                         EMN1      EMN2      EMN3      EMN4     EMN5      IMN1
## CR7d.CTLvsCR7d.CKO 0.3964044 0.4170971 0.5165039 0.6549232 0.469626 0.2860499
##                          IMN2      IMN3      IMN4         IN1       IN2
## CR7d.CTLvsCR7d.CKO 0.01985923 0.5701387 0.9166744 0.001347195 0.7900928
##                          IN3     IPAN1     IPAN2      IPAN3    IPAN4
## CR7d.CTLvsCR7d.CKO 0.4605807 0.1626049 0.6988405 0.09326914 0.285167
#options (scipen = 999)
round(glm.nb_anova.2.intAnno1_df,4)
##                      EMN1   EMN2   EMN3   EMN4   EMN5  IMN1   IMN2   IMN3
## CR7d.CTLvsCR7d.CKO 0.3964 0.4171 0.5165 0.6549 0.4696 0.286 0.0199 0.5701
##                      IMN4    IN1    IN2    IN3  IPAN1  IPAN2  IPAN3  IPAN4
## CR7d.CTLvsCR7d.CKO 0.9167 0.0013 0.7901 0.4606 0.1626 0.6988 0.0933 0.2852
Sort.N <- c("CR7d.CTL","CR7d.CKO")

glm.nb_anova.2.intAnno2 <- list()

for(x1 in 1:2){
  for(x2 in 1:2){
    if(x2>x1){
      
      stat2_intAnno2.s_N <- stat2_intAnno2.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
      
      total.N <- stat2_intAnno2.s_N %>%
        group_by(cnt, rep) %>%
        summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
      
      rownames(total.N) <- paste0(as.character(total.N$cnt),
                                  "_",
                                  as.character(total.N$rep))
      
      stat2_intAnno2.s_N$total <- total.N[paste0(stat2_intAnno2.s_N$cnt,
                                            "_",
                                            stat2_intAnno2.s_N$rep),"total"]
      
      glm.nb_anova.2.intAnno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
      
      for(AA in levels(stat2_intAnno2.s_N$intAnno2)){
        glm.nb_anova.2.intAnno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
          aaa <- stat2_intAnno2.s_N %>% filter(intAnno2==AA)
          aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
          aaa.anova <- anova(aaa.glmnb, test = "LRT")
          aaa.anova$'Pr(>Chi)'[2]
        }, error=function(e){
        
          tryCatch({
            aaa <- stat2_intAnno2.s_N %>% filter(intAnno2==AA)
            aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
            aaa.anova <- anova(aaa.glmnb, test = "LRT")
            aaa.anova$'Pr(>Chi)'[2]
          }, error=function(e){
            
            tryCatch({
              aaa <- stat2_intAnno2.s_N %>% filter(intAnno2==AA)
              aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
              aaa.anova <- anova(aaa.glmnb, test = "LRT")
              aaa.anova$'Pr(>Chi)'[2]
            }, error = function(e){
              NA
            })     
          })
        })
      }
      glm.nb_anova.2.intAnno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.2.intAnno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]]) 
    }
  }
}
glm.nb_anova.2.intAnno2_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.2.intAnno2)))
rownames(glm.nb_anova.2.intAnno2_df) <- names(glm.nb_anova.2.intAnno2)
colnames(glm.nb_anova.2.intAnno2_df) <- gsub("X","C",colnames(glm.nb_anova.2.intAnno2_df))
glm.nb_anova.2.intAnno2_df
##                         EMN1      EMN2      EMN3      EMN4     EMN5      IMN1
## CR7d.CTLvsCR7d.CKO 0.3964044 0.4170971 0.5165039 0.6549232 0.469626 0.2860499
##                          IMN2      IMN3      IMN4         IN1       IN2
## CR7d.CTLvsCR7d.CKO 0.01985923 0.5701387 0.9166744 0.001347195 0.7900928
##                          IN3   IPAN1.1   IPAN1.2   IPAN2.1   IPAN2.2      IPAN3
## CR7d.CTLvsCR7d.CKO 0.4605807 0.1757535 0.3257703 0.9831953 0.5659296 0.09326914
##                       IPAN4
## CR7d.CTLvsCR7d.CKO 0.285167
#options (scipen = 999)
round(glm.nb_anova.2.intAnno2_df,4)
##                      EMN1   EMN2   EMN3   EMN4   EMN5  IMN1   IMN2   IMN3
## CR7d.CTLvsCR7d.CKO 0.3964 0.4171 0.5165 0.6549 0.4696 0.286 0.0199 0.5701
##                      IMN4    IN1    IN2    IN3 IPAN1.1 IPAN1.2 IPAN2.1 IPAN2.2
## CR7d.CTLvsCR7d.CKO 0.9167 0.0013 0.7901 0.4606  0.1758  0.3258  0.9832  0.5659
##                     IPAN3  IPAN4
## CR7d.CTLvsCR7d.CKO 0.0933 0.2852

all markers

intAnno1

# find markers for every cluster compared to all remaining cells
Idents(GEX.seur) <- "intAnno1"

GEX.seur <- PrepSCTFindMarkers(GEX.seur, assay = "SCT")
#GEX.markers.pre <- FindAllMarkers(GEX.seur, only.pos = TRUE, min.pct = 0.01,
#                                  assay = "SCT",
#                                  test.use = "MAST",
#                                  logfc.threshold = 0.1)
GEX.markers.pre <- read.table("Baf53cre_CR.markers.SCT_intAnno1.202405.csv", header = TRUE, sep = ",")
#GEX.markers.pre %>% group_by(cluster) %>% top_n(n = 8, wt = avg_log2FC)
GEX.markers.pre$cluster <- factor(as.character(GEX.markers.pre$cluster),
                          levels = levels(GEX.seur$intAnno1))


markers.pre_t60 <- (GEX.markers.pre %>% group_by(cluster) %>% 
                  filter(pct.1>0.05 & gene %in% grep("Rps|Rpl|mt-|Hsp",gene,invert = T,value = T)) %>%
                   top_n(n = 60, wt = avg_log2FC) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene

markers.pre_t120 <- (GEX.markers.pre %>% group_by(cluster) %>% 
                  filter(pct.1>0.05 & gene %in% grep("Rps|Rpl|mt-|Hsp",gene,invert = T,value = T)) %>%
                   top_n(n = 120, wt = avg_log2FC) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene
ttt = 468
ttt/60
## [1] 7.8
ttt/64
## [1] 7.3125
ttt/65
## [1] 7.2
pp.t60 <- list()

for(i in 1:8){
pp.t60[[i]] <- DotPlot(GEX.seur, features = rev(markers.pre_t60[(60*i-59):(60*i)]))  + coord_flip() + 
           theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
}

pp.t60
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

ttt = 839
ttt/60
## [1] 13.98333
ttt/64
## [1] 13.10938
ttt/65
## [1] 12.90769
pp.t120 <- list()

for(i in 1:14){
pp.t120[[i]] <- DotPlot(GEX.seur, features = rev(markers.pre_t120[(60*i-59):(60*i)]))  + coord_flip() + 
           theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
}

pp.t120
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

## 
## [[13]]

## 
## [[14]]

DEGs

head(GEX.seur@meta.data)
##                        orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCAAGAATACAC-1_2 Stst.CTL_CKO       2282         1290  0.1314636  0.1314636
## AAACCCAAGCAATAGT-1_2 Stst.CTL_CKO       1535         1011  0.1302932  0.5863192
## AAACCCAAGGTGAGCT-1_2 Stst.CTL_CKO       1049          737  0.0000000  0.4766444
## AAACCCACAACGAGGT-1_2 Stst.CTL_CKO       1591         1046  0.6913891  0.4399749
## AAACCCACAAGAGTAT-1_2 Stst.CTL_CKO       2106         1250  0.2849003  0.2849003
## AAACCCACAATCAAGA-1_2 Stst.CTL_CKO       1451          953  0.1378360  0.1378360
##                      FB.info     S.Score    G2M.Score Phase      cnt  rep
## AAACCCAAGAATACAC-1_2   CTL.3 -0.00110357  0.012058898   G2M Stst.CTL rep3
## AAACCCAAGCAATAGT-1_2   CKO.3  0.02907907 -0.015558699     S Stst.CKO rep3
## AAACCCAAGGTGAGCT-1_2   CTL.1 -0.01100811 -0.008015087    G1 Stst.CTL rep1
## AAACCCACAACGAGGT-1_2   CTL.1 -0.02143685  0.005086498   G2M Stst.CTL rep1
## AAACCCACAAGAGTAT-1_2   CKO.3  0.02965845 -0.009057774     S Stst.CKO rep3
## AAACCCACAATCAAGA-1_2   CTL.2 -0.01216686  0.004143546   G2M Stst.CTL rep2
##                         sample tissue nCount_SCT nFeature_SCT seurat_clusters
## AAACCCAAGAATACAC-1_2 Stst.CTL3  Ileum       1865         1287               8
## AAACCCAAGCAATAGT-1_2 Stst.CKO3  Ileum       1525         1010               2
## AAACCCAAGGTGAGCT-1_2 Stst.CTL1  Ileum       1170          735               0
## AAACCCACAACGAGGT-1_2 Stst.CTL1  Ileum       1571         1045               3
## AAACCCACAAGAGTAT-1_2 Stst.CKO3  Ileum       1837         1248              16
## AAACCCACAATCAAGA-1_2 Stst.CTL2  Ileum       1452          953               3
##                      sort_clusters intAnno1 intAnno2
## AAACCCAAGAATACAC-1_2             8     EMN1     EMN1
## AAACCCAAGCAATAGT-1_2             2     IMN1     IMN1
## AAACCCAAGGTGAGCT-1_2             0     EMN1     EMN1
## AAACCCACAACGAGGT-1_2             3     EMN1     EMN1
## AAACCCACAAGAGTAT-1_2            16     IMN3     IMN3
## AAACCCACAATCAAGA-1_2             3     EMN1     EMN1
Idents(GEX.seur) <- "cnt"

#GEX.seur <- PrepSCTFindMarkers(GEX.seur, assay = "SCT")
neur.clusters <- grep("EMN|IMN|IN|IPAN",levels(GEX.seur$intAnno2), value = T)
neur.clusters
##  [1] "EMN1"    "EMN2"    "EMN3"    "EMN4"    "EMN5"    "IMN1"    "IMN2"   
##  [8] "IMN3"    "IMN4"    "IN1"     "IN2"     "IN3"     "IPAN1.1" "IPAN1.2"
## [15] "IPAN2.1" "IPAN2.2" "IPAN3"   "IPAN4"
#
all_new <- neur.clusters
all_new
##  [1] "EMN1"    "EMN2"    "EMN3"    "EMN4"    "EMN5"    "IMN1"    "IMN2"   
##  [8] "IMN3"    "IMN4"    "IN1"     "IN2"     "IN3"     "IPAN1.1" "IPAN1.2"
## [15] "IPAN2.1" "IPAN2.2" "IPAN3"   "IPAN4"
#
list_new <- list()
list_new[["All"]] <- all_new

list_new[["EMNs"]] <- grep("EMN",all_new,value = T)
list_new[["IMNs"]] <- grep("IMN",all_new,value = T)

list_new[["IPAN1"]] <- grep("IPAN1",all_new,value = T)
list_new[["IPAN2"]] <- grep("IPAN2",all_new,value = T)

#list_new[["INs"]] <- grep("IN",all_new,value = T)
#list_new[["IPANs"]] <- grep("IPAN",all_new,value = T)

for(NN in all_new){
  list_new[[NN]] <- NN
}

names_new <- names(list_new)
list_new
## $All
##  [1] "EMN1"    "EMN2"    "EMN3"    "EMN4"    "EMN5"    "IMN1"    "IMN2"   
##  [8] "IMN3"    "IMN4"    "IN1"     "IN2"     "IN3"     "IPAN1.1" "IPAN1.2"
## [15] "IPAN2.1" "IPAN2.2" "IPAN3"   "IPAN4"  
## 
## $EMNs
## [1] "EMN1" "EMN2" "EMN3" "EMN4" "EMN5"
## 
## $IMNs
## [1] "IMN1" "IMN2" "IMN3" "IMN4"
## 
## $IPAN1
## [1] "IPAN1.1" "IPAN1.2"
## 
## $IPAN2
## [1] "IPAN2.1" "IPAN2.2"
## 
## $EMN1
## [1] "EMN1"
## 
## $EMN2
## [1] "EMN2"
## 
## $EMN3
## [1] "EMN3"
## 
## $EMN4
## [1] "EMN4"
## 
## $EMN5
## [1] "EMN5"
## 
## $IMN1
## [1] "IMN1"
## 
## $IMN2
## [1] "IMN2"
## 
## $IMN3
## [1] "IMN3"
## 
## $IMN4
## [1] "IMN4"
## 
## $IN1
## [1] "IN1"
## 
## $IN2
## [1] "IN2"
## 
## $IN3
## [1] "IN3"
## 
## $IPAN1.1
## [1] "IPAN1.1"
## 
## $IPAN1.2
## [1] "IPAN1.2"
## 
## $IPAN2.1
## [1] "IPAN2.1"
## 
## $IPAN2.2
## [1] "IPAN2.2"
## 
## $IPAN3
## [1] "IPAN3"
## 
## $IPAN4
## [1] "IPAN4"
names_new
##  [1] "All"     "EMNs"    "IMNs"    "IPAN1"   "IPAN2"   "EMN1"    "EMN2"   
##  [8] "EMN3"    "EMN4"    "EMN5"    "IMN1"    "IMN2"    "IMN3"    "IMN4"   
## [15] "IN1"     "IN2"     "IN3"     "IPAN1.1" "IPAN1.2" "IPAN2.1" "IPAN2.2"
## [22] "IPAN3"   "IPAN4"

Stst

#test1.seur <- subset(GEX.seur, subset= cnt %in% c("Stst.CTL","Stst.CKO"))
test1.seur
## An object of class Seurat 
## 47749 features across 11409 samples within 3 assays 
## Active assay: SCT (21008 features, 0 variable features)
##  2 other assays present: RNA, integrated
##  3 dimensional reductions calculated: pca, tsne, umap
#head(test1.seur@active.ident)
Idents(test1.seur) <- "cnt"

run

#test1.DEGs_new
# save DEGs' table
df_test1.DEGs_new <- data.frame()
for(nn in names_new){

  df_test1.DEGs_new <- rbind(df_test1.DEGs_new, data.frame(test1.DEGs_new[[nn]],intAnno=nn))
  
  
}

#write.csv(df_test1.DEGs_new, "Baf53cre_CR.DEGs.Stst_CTLvsCKO.intAnno2.csv")
#write.csv(df_test1.DEGs_new %>% filter(intAnno %in% setdiff(names_new,paste0("IPAN",c(1.1,1.2,2.1,2.2)))), "Baf53cre_CR.DEGs.Stst_CTLvsCKO.intAnno1.csv")

stat

cut0: padj 0.05 log2FC 0.1
cut1: padj 0.05 log2FC 0.3
cut2: padj 0.01 log2FC > log2(1.5) (|FC|>1.5)

cut0

# cut1
cut.padj = 0.05
cut.log2FC = 0.1  
  
cut.pct1 = 0.1
stat_test1o.DEGs_new <- df_test1.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(intAnno,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame()
stat_test1o.DEGs_new
##    intAnno  cluster up.DEGs
## 1      All Stst.CTL      15
## 2      All Stst.CKO      17
## 3     EMN1 Stst.CTL       7
## 4     EMN1 Stst.CKO       4
## 5     EMN2 Stst.CTL       1
## 6     EMN3 Stst.CTL       1
## 7     EMN4 Stst.CTL       1
## 8     EMN4 Stst.CKO       1
## 9     EMN5 Stst.CTL       1
## 10    EMNs Stst.CTL      11
## 11    EMNs Stst.CKO       8
## 12    IMN1 Stst.CTL       4
## 13    IMN1 Stst.CKO       2
## 14    IMN2 Stst.CTL       1
## 15    IMN3 Stst.CTL       1
## 16    IMNs Stst.CTL       6
## 17    IMNs Stst.CKO       5
## 18     IN1 Stst.CTL       2
## 19     IN2 Stst.CTL       1
## 20   IPAN1 Stst.CTL       4
## 21   IPAN1 Stst.CKO       2
## 22 IPAN1.1 Stst.CTL       1
## 23 IPAN1.1 Stst.CKO       1
## 24 IPAN1.2 Stst.CTL       2
## 25   IPAN2 Stst.CTL       3
## 26   IPAN2 Stst.CKO       2
## 27 IPAN2.1 Stst.CTL       2
## 28 IPAN2.1 Stst.CKO       2
## 29 IPAN2.2 Stst.CTL       1
## 30   IPAN4 Stst.CTL       1
## 31   IPAN4 Stst.CKO       1
df_test1.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(intAnno,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame() %>%
  ggplot(aes(x=intAnno, y=up.DEGs, color = cluster)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.test1, name="") +
  labs(title=paste0("up.DEGs stat, pct.1>",cut.pct1,", padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Count") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
        title =element_text(size=12, face='bold'))

cut1

# cut1
cut.padj = 0.05
cut.log2FC = 0.3  
  
cut.pct1 = 0.1
stat_test1a.DEGs_new <- df_test1.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(intAnno,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame()
stat_test1a.DEGs_new
##    intAnno  cluster up.DEGs
## 1      All Stst.CTL       1
## 2      All Stst.CKO       2
## 3     EMN1 Stst.CTL       4
## 4     EMN1 Stst.CKO       3
## 5     EMN2 Stst.CTL       1
## 6     EMN3 Stst.CTL       1
## 7     EMN4 Stst.CTL       1
## 8     EMN4 Stst.CKO       1
## 9     EMN5 Stst.CTL       1
## 10    EMNs Stst.CTL       2
## 11    EMNs Stst.CKO       3
## 12    IMN1 Stst.CTL       1
## 13    IMN1 Stst.CKO       1
## 14    IMN2 Stst.CTL       1
## 15    IMN3 Stst.CTL       1
## 16    IMNs Stst.CTL       1
## 17    IMNs Stst.CKO       2
## 18     IN1 Stst.CTL       2
## 19     IN2 Stst.CTL       1
## 20   IPAN1 Stst.CTL       4
## 21   IPAN1 Stst.CKO       2
## 22 IPAN1.1 Stst.CTL       1
## 23 IPAN1.1 Stst.CKO       1
## 24 IPAN1.2 Stst.CTL       2
## 25   IPAN2 Stst.CTL       3
## 26   IPAN2 Stst.CKO       2
## 27 IPAN2.1 Stst.CTL       2
## 28 IPAN2.1 Stst.CKO       2
## 29 IPAN2.2 Stst.CTL       1
## 30   IPAN4 Stst.CTL       1
## 31   IPAN4 Stst.CKO       1
df_test1.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(intAnno,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame() %>%
  ggplot(aes(x=intAnno, y=up.DEGs, color = cluster)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.test1, name="") +
  labs(title=paste0("up.DEGs stat, pct.1>",cut.pct1,", padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Count") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
        title =element_text(size=12, face='bold'))

cut2

# cut2
cut.padj = 0.01
cut.log2FC = log2(1.5)  
  
cut.pct1 = 0.1
stat_test1b.DEGs_new <- df_test1.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(intAnno,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame()
stat_test1b.DEGs_new
##    intAnno  cluster up.DEGs
## 1      All Stst.CTL       1
## 2     EMN1 Stst.CTL       1
## 3     EMN2 Stst.CTL       1
## 4     EMN3 Stst.CTL       1
## 5     EMN4 Stst.CTL       1
## 6     EMN4 Stst.CKO       1
## 7     EMN5 Stst.CTL       1
## 8     EMNs Stst.CTL       1
## 9     IMN1 Stst.CTL       1
## 10    IMN2 Stst.CTL       1
## 11    IMN3 Stst.CTL       1
## 12    IMNs Stst.CTL       1
## 13     IN1 Stst.CTL       2
## 14   IPAN1 Stst.CTL       2
## 15   IPAN1 Stst.CKO       1
## 16 IPAN1.1 Stst.CTL       1
## 17 IPAN1.1 Stst.CKO       1
## 18 IPAN1.2 Stst.CTL       1
## 19   IPAN2 Stst.CTL       2
## 20 IPAN2.1 Stst.CTL       2
## 21   IPAN4 Stst.CTL       1
df_test1.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(intAnno,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame() %>%
  ggplot(aes(x=intAnno, y=up.DEGs, color = cluster)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.test1, name="") +
  labs(title=paste0("up.DEGs stat, pct.1>",cut.pct1,", padj<",cut.padj,", |FC|>",2^cut.log2FC), y = "Count") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
        title =element_text(size=12, face='bold'))

plot

pp_test1.DEGs.new <- lapply(list_new, function(x){NA})
#
test1.DEGs_new.combine <- test1.DEGs_new
test1.DEGs_new.combine <- lapply(test1.DEGs_new.combine, function(x){
  x[x$cluster=="Stst.CTL","avg_log2FC"] <- -x[x$cluster=="Stst.CTL","avg_log2FC"]
  x
})

All

pp_test1.DEGs.new$All <- test1.DEGs_new.combine$All %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Stst.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Stst.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Stst.CTL"=as.vector(color.test1[1]),
                               "Stst.CKO"=as.vector(color.test1[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="All Stst.CTL vs Stst.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$All

EMNs

pp_test1.DEGs.new$EMNs <- test1.DEGs_new.combine$EMNs %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Stst.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Stst.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Stst.CTL"=as.vector(color.test1[1]),
                               "Stst.CKO"=as.vector(color.test1[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="EMNs Stst.CTL vs Stst.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$EMNs

EMN1

pp_test1.DEGs.new$EMN1 <- test1.DEGs_new.combine$EMN1 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Stst.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Stst.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Stst.CTL"=as.vector(color.test1[1]),
                               "Stst.CKO"=as.vector(color.test1[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="EMN1 Stst.CTL vs Stst.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$EMN1

EMN2

pp_test1.DEGs.new$EMN2 <- test1.DEGs_new.combine$EMN2 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Stst.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Stst.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Stst.CTL"=as.vector(color.test1[1]),
                               "Stst.CKO"=as.vector(color.test1[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="EMN2 Stst.CTL vs Stst.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$EMN2

EMN3

pp_test1.DEGs.new$EMN3 <- test1.DEGs_new.combine$EMN3 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Stst.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Stst.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Stst.CTL"=as.vector(color.test1[1]),
                               "Stst.CKO"=as.vector(color.test1[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="EMN3 Stst.CTL vs Stst.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$EMN3

EMN4

pp_test1.DEGs.new$EMN4 <- test1.DEGs_new.combine$EMN4 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Stst.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Stst.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Stst.CTL"=as.vector(color.test1[1]),
                               "Stst.CKO"=as.vector(color.test1[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="EMN4 Stst.CTL vs Stst.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$EMN4

EMN5

pp_test1.DEGs.new$EMN5 <- test1.DEGs_new.combine$EMN5 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Stst.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Stst.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Stst.CTL"=as.vector(color.test1[1]),
                               "Stst.CKO"=as.vector(color.test1[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="EMN5 Stst.CTL vs Stst.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$EMN5

IMNs

pp_test1.DEGs.new$IMNs <- test1.DEGs_new.combine$IMNs %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Stst.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Stst.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Stst.CTL"=as.vector(color.test1[1]),
                               "Stst.CKO"=as.vector(color.test1[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IMNs Stst.CTL vs Stst.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IMNs

IMN1

pp_test1.DEGs.new$IMN1 <- test1.DEGs_new.combine$IMN1 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Stst.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Stst.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Stst.CTL"=as.vector(color.test1[1]),
                               "Stst.CKO"=as.vector(color.test1[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IMN1 Stst.CTL vs Stst.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IMN1

IMN2

pp_test1.DEGs.new$IMN2 <- test1.DEGs_new.combine$IMN2 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Stst.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Stst.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Stst.CTL"=as.vector(color.test1[1]),
                               "Stst.CKO"=as.vector(color.test1[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IMN2 Stst.CTL vs Stst.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IMN2

IMN3

pp_test1.DEGs.new$IMN3 <- test1.DEGs_new.combine$IMN3 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Stst.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Stst.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Stst.CTL"=as.vector(color.test1[1]),
                               "Stst.CKO"=as.vector(color.test1[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IMN3 Stst.CTL vs Stst.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IMN3

IMN4

pp_test1.DEGs.new$IMN4 <- test1.DEGs_new.combine$IMN4 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Stst.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Stst.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Stst.CTL"=as.vector(color.test1[1]),
                               "Stst.CKO"=as.vector(color.test1[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IMN4 Stst.CTL vs Stst.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IMN4

IN1

pp_test1.DEGs.new$IN1 <- test1.DEGs_new.combine$IN1 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Stst.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Stst.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Stst.CTL"=as.vector(color.test1[1]),
                               "Stst.CKO"=as.vector(color.test1[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IN1 Stst.CTL vs Stst.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IN1

IN2

pp_test1.DEGs.new$IN2 <- test1.DEGs_new.combine$IN2 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Stst.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Stst.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Stst.CTL"=as.vector(color.test1[1]),
                               "Stst.CKO"=as.vector(color.test1[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IN2 Stst.CTL vs Stst.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IN2

IN3

pp_test1.DEGs.new$IN3 <- test1.DEGs_new.combine$IN3 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Stst.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Stst.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Stst.CTL"=as.vector(color.test1[1]),
                               "Stst.CKO"=as.vector(color.test1[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IN3 Stst.CTL vs Stst.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IN3

IPAN1

pp_test1.DEGs.new$IPAN1 <- test1.DEGs_new.combine$IPAN1 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Stst.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Stst.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Stst.CTL"=as.vector(color.test1[1]),
                               "Stst.CKO"=as.vector(color.test1[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN1 Stst.CTL vs Stst.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IPAN1

IPAN1.1

pp_test1.DEGs.new$IPAN1.1 <- test1.DEGs_new.combine$IPAN1.1 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Stst.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Stst.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Stst.CTL"=as.vector(color.test1[1]),
                               "Stst.CKO"=as.vector(color.test1[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN1.1 Stst.CTL vs Stst.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IPAN1.1

IPAN1.2

pp_test1.DEGs.new$IPAN1.2 <- test1.DEGs_new.combine$IPAN1.2 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Stst.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Stst.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Stst.CTL"=as.vector(color.test1[1]),
                               "Stst.CKO"=as.vector(color.test1[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN1.2 Stst.CTL vs Stst.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IPAN1.2

IPAN2

pp_test1.DEGs.new$IPAN2 <- test1.DEGs_new.combine$IPAN2 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Stst.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Stst.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Stst.CTL"=as.vector(color.test1[1]),
                               "Stst.CKO"=as.vector(color.test1[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN2 Stst.CTL vs Stst.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IPAN2

IPAN2.1

pp_test1.DEGs.new$IPAN2.1 <- test1.DEGs_new.combine$IPAN2.1 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Stst.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Stst.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Stst.CTL"=as.vector(color.test1[1]),
                               "Stst.CKO"=as.vector(color.test1[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN2.1 Stst.CTL vs Stst.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IPAN2.1

IPAN2.2

pp_test1.DEGs.new$IPAN2.2 <- test1.DEGs_new.combine$IPAN2.2 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Stst.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Stst.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Stst.CTL"=as.vector(color.test1[1]),
                               "Stst.CKO"=as.vector(color.test1[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN2.2 Stst.CTL vs Stst.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IPAN2.2

IPAN3

pp_test1.DEGs.new$IPAN3 <- test1.DEGs_new.combine$IPAN3 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Stst.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Stst.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Stst.CTL"=as.vector(color.test1[1]),
                               "Stst.CKO"=as.vector(color.test1[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN3 Stst.CTL vs Stst.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IPAN3

IPAN4

pp_test1.DEGs.new$IPAN4 <- test1.DEGs_new.combine$IPAN4 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Stst.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Stst.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Stst.CTL"=as.vector(color.test1[1]),
                               "Stst.CKO"=as.vector(color.test1[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN4 Stst.CTL vs Stst.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IPAN4

CR7d

#test2.seur <- subset(GEX.seur, subset= cnt %in% c("CR7d.CTL","CR7d.CKO"))
test2.seur
## An object of class Seurat 
## 47749 features across 12899 samples within 3 assays 
## Active assay: SCT (21008 features, 0 variable features)
##  2 other assays present: RNA, integrated
##  3 dimensional reductions calculated: pca, tsne, umap
#head(test2.seur@active.ident)
Idents(test2.seur) <- "cnt"

run

#test2.DEGs_new
# save DEGs' table
df_test2.DEGs_new <- data.frame()
for(nn in names_new){
  df_test2.DEGs_new <- rbind(df_test2.DEGs_new, data.frame(test2.DEGs_new[[nn]],intAnno=nn))
}

#write.csv(df_test2.DEGs_new, "Baf53cre_CR.DEGs.CR7d_CTLvsCKO.intAnno2.csv")
#write.csv(df_test2.DEGs_new %>% filter(intAnno %in% setdiff(names_new,paste0("IPAN",c(1.1,1.2,2.1,2.2)))), "Baf53cre_CR.DEGs.CR7d_CTLvsCKO.intAnno1.csv")

stat

cut0: padj 0.05 log2FC 0.1
cut1: padj 0.05 log2FC 0.3
cut2: padj 0.01 log2FC > log2(1.5) (|FC|>1.5)

cut0

# cut1
cut.padj = 0.05
cut.log2FC = 0.1  
  
cut.pct1 = 0.1
stat_test2o.DEGs_new <- df_test2.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(intAnno,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame()
stat_test2o.DEGs_new
##    intAnno  cluster up.DEGs
## 1      All CR7d.CTL       4
## 2      All CR7d.CKO      12
## 3     EMN1 CR7d.CTL       3
## 4     EMN1 CR7d.CKO       3
## 5     EMN2 CR7d.CTL       2
## 6     EMN3 CR7d.CTL       2
## 7     EMN4 CR7d.CTL       2
## 8     EMN4 CR7d.CKO       2
## 9     EMN5 CR7d.CTL       2
## 10    EMNs CR7d.CTL       3
## 11    EMNs CR7d.CKO       6
## 12    IMN1 CR7d.CTL       3
## 13    IMN1 CR7d.CKO       4
## 14    IMN2 CR7d.CTL       2
## 15    IMN3 CR7d.CTL       1
## 16    IMNs CR7d.CTL       4
## 17    IMNs CR7d.CKO       4
## 18     IN1 CR7d.CTL       1
## 19     IN1 CR7d.CKO       1
## 20   IPAN1 CR7d.CTL       1
## 21   IPAN1 CR7d.CKO       1
## 22 IPAN1.1 CR7d.CTL       1
## 23 IPAN1.2 CR7d.CTL       1
## 24 IPAN1.2 CR7d.CKO       1
## 25   IPAN2 CR7d.CTL       1
## 26   IPAN2 CR7d.CKO       3
## 27 IPAN2.1 CR7d.CTL       1
## 28 IPAN2.1 CR7d.CKO       1
## 29 IPAN2.2 CR7d.CTL       1
## 30 IPAN2.2 CR7d.CKO       1
df_test2.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(intAnno,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame() %>%
  ggplot(aes(x=intAnno, y=up.DEGs, color = cluster)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.test2, name="") +
  labs(title=paste0("up.DEGs stat, pct.1>",cut.pct1,", padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Count") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
        title =element_text(size=12, face='bold'))

cut1

# cut1
cut.padj = 0.05
cut.log2FC = 0.3  
  
cut.pct1 = 0.1
stat_test2a.DEGs_new <- df_test2.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(intAnno,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame()
stat_test2a.DEGs_new
##    intAnno  cluster up.DEGs
## 1      All CR7d.CTL       1
## 2      All CR7d.CKO       2
## 3     EMN1 CR7d.CTL       3
## 4     EMN1 CR7d.CKO       2
## 5     EMN2 CR7d.CTL       1
## 6     EMN3 CR7d.CTL       2
## 7     EMN4 CR7d.CTL       2
## 8     EMN4 CR7d.CKO       2
## 9     EMN5 CR7d.CTL       2
## 10    EMNs CR7d.CTL       3
## 11    EMNs CR7d.CKO       2
## 12    IMN1 CR7d.CTL       1
## 13    IMN1 CR7d.CKO       4
## 14    IMN2 CR7d.CTL       2
## 15    IMN3 CR7d.CTL       1
## 16    IMNs CR7d.CTL       2
## 17    IMNs CR7d.CKO       4
## 18     IN1 CR7d.CTL       1
## 19     IN1 CR7d.CKO       1
## 20   IPAN1 CR7d.CTL       1
## 21   IPAN1 CR7d.CKO       1
## 22 IPAN1.1 CR7d.CTL       1
## 23 IPAN1.2 CR7d.CTL       1
## 24 IPAN1.2 CR7d.CKO       1
## 25   IPAN2 CR7d.CTL       1
## 26   IPAN2 CR7d.CKO       3
## 27 IPAN2.1 CR7d.CTL       1
## 28 IPAN2.1 CR7d.CKO       1
## 29 IPAN2.2 CR7d.CTL       1
## 30 IPAN2.2 CR7d.CKO       1
df_test2.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(intAnno,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame() %>%
  ggplot(aes(x=intAnno, y=up.DEGs, color = cluster)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.test2, name="") +
  labs(title=paste0("up.DEGs stat, pct.1>",cut.pct1,", padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Count") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
        title =element_text(size=12, face='bold'))

cut2

# cut2
cut.padj = 0.01
cut.log2FC = log2(1.5)  
  
cut.pct1 = 0.1
stat_test2b.DEGs_new <- df_test2.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(intAnno,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame()
stat_test2b.DEGs_new
##    intAnno  cluster up.DEGs
## 1      All CR7d.CTL       1
## 2     EMN1 CR7d.CTL       1
## 3     EMN2 CR7d.CTL       1
## 4     EMN3 CR7d.CTL       1
## 5     EMN4 CR7d.CTL       1
## 6     EMN5 CR7d.CTL       2
## 7     EMNs CR7d.CTL       1
## 8     IMN1 CR7d.CTL       1
## 9     IMN2 CR7d.CTL       1
## 10    IMNs CR7d.CTL       1
## 11     IN1 CR7d.CTL       1
## 12     IN1 CR7d.CKO       1
## 13   IPAN1 CR7d.CTL       1
## 14   IPAN1 CR7d.CKO       1
## 15 IPAN1.1 CR7d.CTL       1
## 16 IPAN1.2 CR7d.CTL       1
## 17 IPAN1.2 CR7d.CKO       1
## 18   IPAN2 CR7d.CTL       1
## 19   IPAN2 CR7d.CKO       2
## 20 IPAN2.1 CR7d.CTL       1
## 21 IPAN2.2 CR7d.CTL       1
## 22 IPAN2.2 CR7d.CKO       1
df_test2.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(intAnno,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame() %>%
  ggplot(aes(x=intAnno, y=up.DEGs, color = cluster)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.test2, name="") +
  labs(title=paste0("up.DEGs stat, pct.1>",cut.pct1,", padj<",cut.padj,", |FC|>",2^cut.log2FC), y = "Count") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
        title =element_text(size=12, face='bold'))

plot

pp_test2.DEGs.new <- lapply(list_new, function(x){NA})
#
test2.DEGs_new.combine <- test2.DEGs_new
test2.DEGs_new.combine <- lapply(test2.DEGs_new.combine, function(x){
  x[x$cluster=="CR7d.CTL","avg_log2FC"] <- -x[x$cluster=="CR7d.CTL","avg_log2FC"]
  x
})

All

pp_test2.DEGs.new$All <- test2.DEGs_new.combine$All %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CR7d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CR7d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CR7d.CTL"=as.vector(color.test2[1]),
                               "CR7d.CKO"=as.vector(color.test2[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="All CR7d.CTL vs CR7d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$All

EMNs

pp_test2.DEGs.new$EMNs <- test2.DEGs_new.combine$EMNs %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CR7d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CR7d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CR7d.CTL"=as.vector(color.test2[1]),
                               "CR7d.CKO"=as.vector(color.test2[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="EMNs CR7d.CTL vs CR7d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$EMNs

EMN1

pp_test2.DEGs.new$EMN1 <- test2.DEGs_new.combine$EMN1 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CR7d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CR7d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CR7d.CTL"=as.vector(color.test2[1]),
                               "CR7d.CKO"=as.vector(color.test2[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="EMN1 CR7d.CTL vs CR7d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$EMN1

EMN2

pp_test2.DEGs.new$EMN2 <- test2.DEGs_new.combine$EMN2 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CR7d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CR7d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CR7d.CTL"=as.vector(color.test2[1]),
                               "CR7d.CKO"=as.vector(color.test2[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="EMN2 CR7d.CTL vs CR7d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$EMN2

EMN3

pp_test2.DEGs.new$EMN3 <- test2.DEGs_new.combine$EMN3 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CR7d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CR7d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CR7d.CTL"=as.vector(color.test2[1]),
                               "CR7d.CKO"=as.vector(color.test2[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="EMN3 CR7d.CTL vs CR7d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$EMN3

EMN4

pp_test2.DEGs.new$EMN4 <- test2.DEGs_new.combine$EMN4 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CR7d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CR7d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CR7d.CTL"=as.vector(color.test2[1]),
                               "CR7d.CKO"=as.vector(color.test2[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="EMN4 CR7d.CTL vs CR7d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$EMN4

EMN5

pp_test2.DEGs.new$EMN5 <- test2.DEGs_new.combine$EMN5 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CR7d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CR7d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CR7d.CTL"=as.vector(color.test2[1]),
                               "CR7d.CKO"=as.vector(color.test2[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="EMN5 CR7d.CTL vs CR7d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$EMN5

IMNs

pp_test2.DEGs.new$IMNs <- test2.DEGs_new.combine$IMNs %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CR7d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CR7d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CR7d.CTL"=as.vector(color.test2[1]),
                               "CR7d.CKO"=as.vector(color.test2[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IMNs CR7d.CTL vs CR7d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IMNs

IMN1

pp_test2.DEGs.new$IMN1 <- test2.DEGs_new.combine$IMN1 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CR7d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CR7d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CR7d.CTL"=as.vector(color.test2[1]),
                               "CR7d.CKO"=as.vector(color.test2[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IMN1 CR7d.CTL vs CR7d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IMN1

IMN2

pp_test2.DEGs.new$IMN2 <- test2.DEGs_new.combine$IMN2 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CR7d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CR7d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CR7d.CTL"=as.vector(color.test2[1]),
                               "CR7d.CKO"=as.vector(color.test2[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IMN2 CR7d.CTL vs CR7d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IMN2

IMN3

pp_test2.DEGs.new$IMN3 <- test2.DEGs_new.combine$IMN3 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CR7d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CR7d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CR7d.CTL"=as.vector(color.test2[1]),
                               "CR7d.CKO"=as.vector(color.test2[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IMN3 CR7d.CTL vs CR7d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IMN3

IMN4

pp_test2.DEGs.new$IMN4 <- test2.DEGs_new.combine$IMN4 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CR7d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CR7d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CR7d.CTL"=as.vector(color.test2[1]),
                               "CR7d.CKO"=as.vector(color.test2[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IMN4 CR7d.CTL vs CR7d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IMN4

IN1

pp_test2.DEGs.new$IN1 <- test2.DEGs_new.combine$IN1 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CR7d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CR7d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CR7d.CTL"=as.vector(color.test2[1]),
                               "CR7d.CKO"=as.vector(color.test2[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IN1 CR7d.CTL vs CR7d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IN1

IN2

pp_test2.DEGs.new$IN2 <- test2.DEGs_new.combine$IN2 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CR7d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CR7d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CR7d.CTL"=as.vector(color.test2[1]),
                               "CR7d.CKO"=as.vector(color.test2[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IN2 CR7d.CTL vs CR7d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IN2

IN3

pp_test2.DEGs.new$IN3 <- test2.DEGs_new.combine$IN3 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CR7d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CR7d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CR7d.CTL"=as.vector(color.test2[1]),
                               "CR7d.CKO"=as.vector(color.test2[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IN3 CR7d.CTL vs CR7d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IN3

IPAN1

pp_test2.DEGs.new$IPAN1 <- test2.DEGs_new.combine$IPAN1 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CR7d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CR7d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CR7d.CTL"=as.vector(color.test2[1]),
                               "CR7d.CKO"=as.vector(color.test2[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN1 CR7d.CTL vs CR7d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IPAN1

IPAN1.1

pp_test2.DEGs.new$IPAN1.1 <- test2.DEGs_new.combine$IPAN1.1 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CR7d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CR7d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CR7d.CTL"=as.vector(color.test2[1]),
                               "CR7d.CKO"=as.vector(color.test2[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN1.1 CR7d.CTL vs CR7d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IPAN1.1

IPAN1.2

pp_test2.DEGs.new$IPAN1.2 <- test2.DEGs_new.combine$IPAN1.2 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CR7d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CR7d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CR7d.CTL"=as.vector(color.test2[1]),
                               "CR7d.CKO"=as.vector(color.test2[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN1.2 CR7d.CTL vs CR7d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IPAN1.2

IPAN2

pp_test2.DEGs.new$IPAN2 <- test2.DEGs_new.combine$IPAN2 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CR7d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CR7d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CR7d.CTL"=as.vector(color.test2[1]),
                               "CR7d.CKO"=as.vector(color.test2[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN2 CR7d.CTL vs CR7d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IPAN2

IPAN2.1

pp_test2.DEGs.new$IPAN2.1 <- test2.DEGs_new.combine$IPAN2.1 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CR7d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CR7d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CR7d.CTL"=as.vector(color.test2[1]),
                               "CR7d.CKO"=as.vector(color.test2[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN2.1 CR7d.CTL vs CR7d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IPAN2.1

IPAN2.2

pp_test2.DEGs.new$IPAN2.2 <- test2.DEGs_new.combine$IPAN2.2 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CR7d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CR7d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CR7d.CTL"=as.vector(color.test2[1]),
                               "CR7d.CKO"=as.vector(color.test2[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN2.2 CR7d.CTL vs CR7d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IPAN2.2

IPAN3

pp_test2.DEGs.new$IPAN3 <- test2.DEGs_new.combine$IPAN3 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CR7d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CR7d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CR7d.CTL"=as.vector(color.test2[1]),
                               "CR7d.CKO"=as.vector(color.test2[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN3 CR7d.CTL vs CR7d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IPAN3

IPAN4

pp_test2.DEGs.new$IPAN4 <- test2.DEGs_new.combine$IPAN4 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CR7d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CR7d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CR7d.CTL"=as.vector(color.test2[1]),
                               "CR7d.CKO"=as.vector(color.test2[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN4 CR7d.CTL vs CR7d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IPAN4

final markers

short_old

markers.old.s <- list(EMN=c("Chat","Bnc2",#"Tox","Ptprt",
                       "Gfra2","Oprk1",#"Adamtsl1", 
                       "Fbxw15","Fbxw24",#"Chrna7",
                       "Satb1","Cntnap5b",
                       "Gabrb1","Nxph1","Lama2","Lrrc7",
                       "Ryr3",#"Eda",
                       "Tac1",
                       #"Kctd8","Ntrk2",
                       "Penk", 
                       "Fut9","Nfatc1","Egfr","Ahr"#,#"Mgll",
                       #"Chrm3"
                       ),
                 IMN=c("Nos1","Kcnab1",
                       "Gfra1","Etv1",
                       #"Man1a","Airn",
                       "Adcy2","Cmah","Creb5","Vip","Pde1a",
                       "Ebf1"#,"Gpc5"
                       ),
                 IN=c("Npas3","Synpr","St18","Gal",
                      "Neurod6",
                      #"Kcnk13",
                      "Moxd1","Sctr",
                      "Piezo1","Sst",#"Adamts9",
                      "Kcnn2"),
                 IPAN=c("Calb2","Calcb","Nmu","Adgrg6",#"Pcdh10",
                        "Ngfr","Galr1","Il7",#"Aff2",
                        #"Gpr149",
                        "Cdh6","Cdh8",
                        "Clstn2",#"Ano2","Ntrk3",
                        "Cpne4",#"Vwc2l",
                        "Cdh9","Scgn",
                        #"Vcan",
                        "Cck","Piezo2","Kcnh7",
                        #"Rerg",
                        "Bmpr1b","Skap1","Ntng1",
                        "Tafa2","Nxph2"))
pn.intAnno1.test0a <- DotPlot(GEX.seur, features = as.vector(unlist(markers.old.s)), group.by = "intAnno1",
        cols = c("midnightblue","darkorange1")) +
  #coord_flip() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + 
  scale_y_discrete(limits=rev) + 
  labs(title="All intAnno1")
pn.intAnno1.test0a

pn.intAnno2.test0a <- DotPlot(GEX.seur, features = as.vector(unlist(markers.old.s)), group.by = "intAnno2",
        cols = c("midnightblue","darkorange1")) +
  #coord_flip() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + 
  scale_y_discrete(limits=rev) + 
  labs(title="All intAnno2")
pn.intAnno2.test0a

pn.intAnno1.test0b <- DotPlot(GEX.seur, features = as.vector(unlist(markers.old.s)), group.by = "intAnno1",
        cols = c("midnightblue","darkorange1")) +
  coord_flip() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.4,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + scale_y_discrete(limits=rev) + 
  labs(title="All intAnno1")
pn.intAnno1.test0b

pn.intAnno2.test0b <- DotPlot(GEX.seur, features = as.vector(unlist(markers.old.s)), group.by = "intAnno2",
        cols = c("midnightblue","darkorange1")) +
  coord_flip() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.4,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + scale_y_discrete(limits=rev) + 
  labs(title="All intAnno2")
pn.intAnno2.test0b

pn.intAnno1.test1a <- DotPlot(test1.seur, features = as.vector(unlist(markers.old.s)), group.by = "intAnno1",
        cols = c("midnightblue","darkorange1")) +
  #coord_flip() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + 
  scale_y_discrete(limits=rev) + 
  labs(title="Stst intAnno1")
pn.intAnno1.test1a

pn.intAnno2.test1a <- DotPlot(test1.seur, features = as.vector(unlist(markers.old.s)), group.by = "intAnno2",
        cols = c("midnightblue","darkorange1")) +
  #coord_flip() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + 
  scale_y_discrete(limits=rev) + 
  labs(title="Stst intAnno2")
pn.intAnno2.test1a

pn.intAnno1.test1b <- DotPlot(test1.seur, features = as.vector(unlist(markers.old.s)), group.by = "intAnno1",
        cols = c("midnightblue","darkorange1")) +
  coord_flip() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.4,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + scale_y_discrete(limits=rev) + 
  labs(title="Stst intAnno1")
pn.intAnno1.test1b

pn.intAnno2.test1b <- DotPlot(test1.seur, features = as.vector(unlist(markers.old.s)), group.by = "intAnno2",
        cols = c("midnightblue","darkorange1")) +
  coord_flip() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.4,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + scale_y_discrete(limits=rev) + 
  labs(title="Stst intAnno2")
pn.intAnno2.test1b

pn.intAnno1.test2a <- DotPlot(test2.seur, features = as.vector(unlist(markers.old.s)), group.by = "intAnno1",
        cols = c("midnightblue","darkorange1")) +
  #coord_flip() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + 
  scale_y_discrete(limits=rev) + 
  labs(title="CR7d intAnno1")
pn.intAnno1.test2a

pn.intAnno2.test2a <- DotPlot(test2.seur, features = as.vector(unlist(markers.old.s)), group.by = "intAnno2",
        cols = c("midnightblue","darkorange1")) +
  #coord_flip() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + 
  scale_y_discrete(limits=rev) + 
  labs(title="CR7d intAnno2")
pn.intAnno2.test2a

pn.intAnno1.test2b <- DotPlot(test2.seur, features = as.vector(unlist(markers.old.s)), group.by = "intAnno1",
        cols = c("midnightblue","darkorange1")) +
  coord_flip() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.4,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + scale_y_discrete(limits=rev) + 
  labs(title="CR7d intAnno1")
pn.intAnno1.test2b

pn.intAnno2.test2b <- DotPlot(test2.seur, features = as.vector(unlist(markers.old.s)), group.by = "intAnno2",
        cols = c("midnightblue","darkorange1")) +
  coord_flip() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.4,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + scale_y_discrete(limits=rev) + 
  labs(title="CR7d intAnno2")
pn.intAnno2.test2b

ref.seur <- readRDS("../../20230704_10x_SZJ/analysis_ref/GSE149524.P21.integration_Anno.s.rds")
ref.seur
## An object of class Seurat 
## 37583 features across 4419 samples within 3 assays 
## Active assay: SCT (16365 features, 0 variable features)
##  2 other assays present: RNA, integrated
##  3 dimensional reductions calculated: pca, tsne, umap
DimPlot(ref.seur, reduction = "umap", label = T, group.by = "Anno1", cols = color.ref) +
  DimPlot(ref.seur, reduction = "umap", label = T, group.by = "Anno2")

pn.ref.a <- DotPlot(ref.seur, features = as.vector(unlist(markers.old.s)), group.by = "Anno2",
        cols = c("midnightblue","darkorange1")) +
  #coord_flip() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + 
  scale_y_discrete(limits=rev) + 
  labs(title="NatNeur2021 P21")
pn.ref.a

pn.ref.b <- DotPlot(ref.seur, features = as.vector(unlist(markers.old.s)), group.by = "Anno2",
        cols = c("midnightblue","darkorange1")) +
  coord_flip() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.4,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + scale_y_discrete(limits=rev) + 
  labs(title="NatNeur2021 P21")
pn.ref.b

short_new

markers.new.s <- list(EMN=c("Chat","Bnc2",#"Tox","Ptprt",
                       "Gfra2","Oprk1",#"Adamtsl1", 
                       "Fbxw15","Fbxw24",#"Chrna7",
                       "Satb1","Itga6","Cntnap5b",
                       "Chgb","Nxph1",
                       "Lama2","Efnb2","Itgb8",
                       "Lrrc7",
                       "Ryr3",#"Eda",
                       "Tac1",
                       #"Kctd8","Ntrk2",
                       "Penk",
                       "Fut9","Nfatc1","Egfr","Ahr"#"Mgll",
                       #"Chrm3"
                       ),
                 IMN=c("Nos1","Kcnab1",
                       "Gfra1","Etv1",
                       #"Man1a","Airn",
                       "Adcy2","Cmah","Col25a1",
                       "Mid1","Creb5","Vip","Pde1a",
                       "Ebf1",#,"Gpc5"
                       "Ppara","Pcdh11x",
                       "Adcy8","Grp"
                       ),
                 IN=c("Npas3","Synpr","St18","Gal",
                      "Cdh10","Neurod6",
                      "Kcnk13",
                      "Moxd1","Sctr",
                      "Piezo1","Vipr2","Sst",#"Adamts9",
                      "Kcnn2"
                      ),
                 IPAN=c("Calb2","Adcy1",
                        "Nmu","Adgrg6",#"Pcdh10",
                        "Ngfr","Il7",
                        "Itgb6","Calcb","Galr1",
                        #"Aff2",
                        #"Gpr149",
                        "Met",
                        "Cpne4","Cdh6","Cdh8",
                        "Clstn2",#"Ano2","Ntrk3",
                        #"Vwc2l",
                        "Car10","Scgn","Glp2r","Cck",
                        "Cdh9",
                        #"Vcan",
                        "Dcc",
                        "Gabrb1",
                        "Piezo2","Kcnh7",
                        #"Rerg",
                        "Bmpr1b","Ntng1","Skap1",
                        "Tafa2","Nxph2"))
pm.intAnno1.test0a <- DotPlot(GEX.seur, features = as.vector(unlist(markers.new.s)), group.by = "intAnno1",
        cols = c("midnightblue","darkorange1")) +
  #coord_flip() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + 
  scale_y_discrete(limits=rev) + 
  labs(title="All intAnno1")
pm.intAnno1.test0a

pm.intAnno2.test0a <- DotPlot(GEX.seur, features = as.vector(unlist(markers.new.s)), group.by = "intAnno2",
        cols = c("midnightblue","darkorange1")) +
  #coord_flip() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + 
  scale_y_discrete(limits=rev) + 
  labs(title="All intAnno2")
pm.intAnno2.test0a

pm.intAnno1.test0b <- DotPlot(GEX.seur, features = as.vector(unlist(markers.new.s)), group.by = "intAnno1",
        cols = c("midnightblue","darkorange1")) +
  coord_flip() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.4,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + scale_y_discrete(limits=rev) + 
  labs(title="All intAnno1")
pm.intAnno1.test0b

pm.intAnno2.test0b <- DotPlot(GEX.seur, features = as.vector(unlist(markers.new.s)), group.by = "intAnno2",
        cols = c("midnightblue","darkorange1")) +
  coord_flip() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.4,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + scale_y_discrete(limits=rev) + 
  labs(title="All intAnno2")
pm.intAnno2.test0b

pm.intAnno1.test1a <- DotPlot(test1.seur, features = as.vector(unlist(markers.new.s)), group.by = "intAnno1",
        cols = c("midnightblue","darkorange1")) +
  #coord_flip() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + 
  scale_y_discrete(limits=rev) + 
  labs(title="Stst intAnno1")
pm.intAnno1.test1a

pm.intAnno2.test1a <- DotPlot(test1.seur, features = as.vector(unlist(markers.new.s)), group.by = "intAnno2",
        cols = c("midnightblue","darkorange1")) +
  #coord_flip() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + 
  scale_y_discrete(limits=rev) + 
  labs(title="Stst intAnno2")
pm.intAnno2.test1a

pm.intAnno1.test1b <- DotPlot(test1.seur, features = as.vector(unlist(markers.new.s)), group.by = "intAnno1",
        cols = c("midnightblue","darkorange1")) +
  coord_flip() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.4,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + scale_y_discrete(limits=rev) + 
  labs(title="Stst intAnno1")
pm.intAnno1.test1b

pm.intAnno2.test1b <- DotPlot(test1.seur, features = as.vector(unlist(markers.new.s)), group.by = "intAnno2",
        cols = c("midnightblue","darkorange1")) +
  coord_flip() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.4,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + scale_y_discrete(limits=rev) + 
  labs(title="Stst intAnno2")
pm.intAnno2.test1b

pm.intAnno1.test2a <- DotPlot(test2.seur, features = as.vector(unlist(markers.new.s)), group.by = "intAnno1",
        cols = c("midnightblue","darkorange1")) +
  #coord_flip() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + 
  scale_y_discrete(limits=rev) + 
  labs(title="CR7d intAnno1")
pm.intAnno1.test2a

pm.intAnno2.test2a <- DotPlot(test2.seur, features = as.vector(unlist(markers.new.s)), group.by = "intAnno2",
        cols = c("midnightblue","darkorange1")) +
  #coord_flip() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + 
  scale_y_discrete(limits=rev) + 
  labs(title="CR7d intAnno2")
pm.intAnno2.test2a

pm.intAnno1.test2b <- DotPlot(test2.seur, features = as.vector(unlist(markers.new.s)), group.by = "intAnno1",
        cols = c("midnightblue","darkorange1")) +
  coord_flip() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.4,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + scale_y_discrete(limits=rev) + 
  labs(title="CR7d intAnno1")
pm.intAnno1.test2b

pm.intAnno2.test2b <- DotPlot(test2.seur, features = as.vector(unlist(markers.new.s)), group.by = "intAnno2",
        cols = c("midnightblue","darkorange1")) +
  coord_flip() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.4,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + scale_y_discrete(limits=rev) + 
  labs(title="CR7d intAnno2")
pm.intAnno2.test2b

pm.ref.a <- DotPlot(ref.seur, features = as.vector(unlist(markers.new.s)), group.by = "Anno2",
        cols = c("midnightblue","darkorange1")) +
  #coord_flip() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + 
  scale_y_discrete(limits=rev) + 
  labs(title="NatNeur2021 P21")
pm.ref.a

pm.ref.b <- DotPlot(ref.seur, features = as.vector(unlist(markers.new.s)), group.by = "Anno2",
        cols = c("midnightblue","darkorange1")) +
  coord_flip() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.4,size = 9.15))+ 
  scale_color_gradientn(colours  = material.heat(100)) +
  #scale_x_discrete(limits=rev) + scale_y_discrete(limits=rev) + 
  labs(title="NatNeur2021 P21")
pm.ref.b

signature scores

top list

use top120-built top markers

load("I:/Shared_win/projects/20231220_10x_SZJ/analysis_plus/top.markerlist.RData")
data.frame(row.names = "topmarkers",lapply(top.list,length))
##            EMN1 EMN2 EMN3 EMN4 EMN5 IMN1 IMN2 IMN3 IMN4 IN1 IN2 IN3 IPAN1.1
## topmarkers   96   32   37   34   70   71   33   49   53  60  56  63      70
##            IPAN1.2 IPAN2.1 IPAN2.2 IPAN3 IPAN4
## topmarkers      51      43      52    52    50
for(nn in neur.clusters){
  
  GEX.seur <- add_geneset_score(obj = GEX.seur,
                                Assay = "SCT",
                                geneset = top.list[[nn]],
                                setname = paste0("score.",nn))
  
}
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data
mapal <- colorRampPalette(RColorBrewer::brewer.pal(4,"Spectral"))(120)
FeaturePlot(GEX.seur, features = paste0("score.",neur.clusters), ncol = 5) &
  scale_color_gradientn(colors = rev(mapal))

mapal <- colorRampPalette(RColorBrewer::brewer.pal(4,"Spectral"))(120)
FeaturePlot(GEX.seur, features = paste0("score.",neur.clusters), ncol = 5, raster = T, pt.size = 3.5) &
  scale_color_gradientn(colors = rev(mapal))

check some markers

Nmu

check.m <- "Nmu"

#
VlnPlot(test1.seur, features = check.m, group.by = "intAnno1", assay = "RNA") + NoLegend() + labs(x="Stst") +
  scale_fill_manual(values = color.A1)

VlnPlot(test1.seur, features = check.m, group.by = "intAnno1", assay = "RNA", split.by = "cnt") + NoLegend() + labs(x="Stst") +
  scale_fill_manual(values = color.test1)

#
VlnPlot(test2.seur, features = check.m, group.by = "intAnno1", assay = "RNA") + NoLegend() + labs(x="CR7d") +
  scale_fill_manual(values = color.A1)

VlnPlot(test2.seur, features = check.m, group.by = "intAnno1", assay = "RNA", split.by = "cnt") + NoLegend() + labs(x="CR7d") +
  scale_fill_manual(values = color.test2)

##
# use coord_cartesian(ylim()) instead of ylim()

VlnPlot(test1.seur, features = check.m, group.by = "cnt", cols = color.test1, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("All - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,4.5)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Stst.CTL","Stst.CKO")),
                               label.y = c(3.5),
                               size=3.5
                               )

#
VlnPlot(subset(test1.seur, subset=intAnno1 %in% c("IPAN1")), features = check.m, group.by = "cnt", cols = color.test1, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IPAN1 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,4.5)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Stst.CTL","Stst.CKO")),
                               label.y = c(3.5),
                               size=3.5
                               ) 

VlnPlot(subset(test1.seur, subset=intAnno2 %in% c("IPAN1.1")), features = check.m, group.by = "cnt", cols = color.test1, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IPAN1.1 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,4.5)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Stst.CTL","Stst.CKO")),
                               label.y = c(3.5),
                               size=3.5
                               )

VlnPlot(subset(test1.seur, subset=intAnno2 %in% c("IPAN1.2")), features = check.m, group.by = "cnt", cols = color.test1, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IPAN1.2 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,4.5)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Stst.CTL","Stst.CKO")),
                               label.y = c(3.5),
                               size=3.5
                               )

##
# use coord_cartesian(ylim()) instead of ylim()

VlnPlot(test2.seur, features = check.m, group.by = "cnt", cols = color.test2, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("All - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,4.5)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("CR7d.CTL","CR7d.CKO")),
                               label.y = c(3.5),
                               size=3.5
                               )

#
VlnPlot(subset(test2.seur, subset=intAnno1 %in% c("IPAN1")), features = check.m, group.by = "cnt", cols = color.test2, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IPAN1 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,4.5)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("CR7d.CTL","CR7d.CKO")),
                               label.y = c(3.5),
                               size=3.5
                               )

VlnPlot(subset(test2.seur, subset=intAnno2 %in% c("IPAN1.1")), features = check.m, group.by = "cnt", cols = color.test2, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IPAN1.1 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,4.5)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("CR7d.CTL","CR7d.CKO")),
                               label.y = c(3.5),
                               size=3.5
                               )

VlnPlot(subset(test2.seur, subset=intAnno2 %in% c("IPAN1.2")), features = check.m, group.by = "cnt", cols = color.test2, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IPAN1.2 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,4.5)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("CR7d.CTL","CR7d.CKO")),
                               label.y = c(3.5),
                               size=3.5
                               )

Calcb

check.m <- "Calcb"

#
cowplot::plot_grid(
VlnPlot(test1.seur, features = check.m, group.by = "intAnno1", assay = "RNA") + NoLegend() + labs(x="Stst") +
  scale_fill_manual(values = color.A1),

VlnPlot(test2.seur, features = check.m, group.by = "intAnno1", assay = "RNA") + NoLegend() + labs(x="CR7d") +
  scale_fill_manual(values = color.A1),

VlnPlot(test1.seur, features = check.m, group.by = "intAnno1", assay = "RNA", split.by = "cnt") + NoLegend() + labs(x="Stst") +
  scale_fill_manual(values = color.test1),

VlnPlot(test2.seur, features = check.m, group.by = "intAnno1", assay = "RNA", split.by = "cnt") + NoLegend() + labs(x="CR7d") +
  scale_fill_manual(values = color.test2),
#
ncol = 2)

##
# use coord_cartesian(ylim()) instead of ylim()

cowplot::plot_grid(
VlnPlot(test1.seur, features = check.m, group.by = "cnt", cols = color.test1, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("All - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,3.5)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Stst.CTL","Stst.CKO")),
                               label.y = c(2.75),
                               size=3.5
                               ),

#
VlnPlot(subset(test1.seur, subset=intAnno1 %in% c("IPAN1")), features = check.m, group.by = "cnt", cols = color.test1, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IPAN1 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,3.5)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Stst.CTL","Stst.CKO")),
                               label.y = c(2.75),
                               size=3.5
                               ) ,

VlnPlot(subset(test1.seur, subset=intAnno2 %in% c("IPAN1.1")), features = check.m, group.by = "cnt", cols = color.test1, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IPAN1.1 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,3.5)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Stst.CTL","Stst.CKO")),
                               label.y = c(2.75),
                               size=3.5
                               ),

VlnPlot(subset(test1.seur, subset=intAnno2 %in% c("IPAN1.2")), features = check.m, group.by = "cnt", cols = color.test1, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IPAN1.2 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,3.5)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Stst.CTL","Stst.CKO")),
                               label.y = c(2.75),
                               size=3.5
                               ),
ncol = 2)

##
# use coord_cartesian(ylim()) instead of ylim()
cowplot::plot_grid(
VlnPlot(test2.seur, features = check.m, group.by = "cnt", cols = color.test2, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("All - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,3.5)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("CR7d.CTL","CR7d.CKO")),
                               label.y = c(2.75),
                               size=3.5
                               ),

#
VlnPlot(subset(test2.seur, subset=intAnno1 %in% c("IPAN1")), features = check.m, group.by = "cnt", cols = color.test2, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IPAN1 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,3.5)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("CR7d.CTL","CR7d.CKO")),
                               label.y = c(2.75),
                               size=3.5
                               ),

VlnPlot(subset(test2.seur, subset=intAnno2 %in% c("IPAN1.1")), features = check.m, group.by = "cnt", cols = color.test2, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IPAN1.1 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,3.5)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("CR7d.CTL","CR7d.CKO")),
                               label.y = c(2.75),
                               size=3.5
                               ),

VlnPlot(subset(test2.seur, subset=intAnno2 %in% c("IPAN1.2")), features = check.m, group.by = "cnt", cols = color.test2, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IPAN1.2 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,3.5)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("CR7d.CTL","CR7d.CKO")),
                               label.y = c(2.75),
                               size=3.5
                               ),
ncol = 2)

##
# use coord_cartesian(ylim()) instead of ylim()

cowplot::plot_grid(
VlnPlot(test1.seur, features = check.m, group.by = "cnt", cols = color.test1, pt.size = 0, assay = "SCT") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("All - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,2)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Stst.CTL","Stst.CKO")),
                               label.y = c(1.5),
                               size=3.5
                               ),

#
VlnPlot(subset(test1.seur, subset=intAnno1 %in% c("IPAN1")), features = check.m, group.by = "cnt", cols = color.test1, pt.size = 0, assay = "SCT") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IPAN1 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,2)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Stst.CTL","Stst.CKO")),
                               label.y = c(1.5),
                               size=3.5
                               ) ,

VlnPlot(subset(test1.seur, subset=intAnno2 %in% c("IPAN1.1")), features = check.m, group.by = "cnt", cols = color.test1, pt.size = 0, assay = "SCT") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IPAN1.1 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,2)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Stst.CTL","Stst.CKO")),
                               label.y = c(1.5),
                               size=3.5
                               ),

VlnPlot(subset(test1.seur, subset=intAnno2 %in% c("IPAN1.2")), features = check.m, group.by = "cnt", cols = color.test1, pt.size = 0, assay = "SCT") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IPAN1.2 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,2)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Stst.CTL","Stst.CKO")),
                               label.y = c(1.5),
                               size=3.5
                               ),
ncol = 2)

##
# use coord_cartesian(ylim()) instead of ylim()
cowplot::plot_grid(
VlnPlot(test2.seur, features = check.m, group.by = "cnt", cols = color.test2, pt.size = 0, assay = "SCT") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("All - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,2)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("CR7d.CTL","CR7d.CKO")),
                               label.y = c(1.45),
                               size=3.5
                               ),

#
VlnPlot(subset(test2.seur, subset=intAnno1 %in% c("IPAN1")), features = check.m, group.by = "cnt", cols = color.test2, pt.size = 0, assay = "SCT") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IPAN1 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,2)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("CR7d.CTL","CR7d.CKO")),
                               label.y = c(1.45),
                               size=3.5
                               ),

VlnPlot(subset(test2.seur, subset=intAnno2 %in% c("IPAN1.1")), features = check.m, group.by = "cnt", cols = color.test2, pt.size = 0, assay = "SCT") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IPAN1.1 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,2)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("CR7d.CTL","CR7d.CKO")),
                               label.y = c(1.25),
                               size=3.5
                               ),

VlnPlot(subset(test2.seur, subset=intAnno2 %in% c("IPAN1.2")), features = check.m, group.by = "cnt", cols = color.test2, pt.size = 0, assay = "SCT") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IPAN1.2 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,2)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("CR7d.CTL","CR7d.CKO")),
                               label.y = c(1.45),
                               size=3.5
                               ),
ncol = 2)

Gal

check.m <- "Gal"

#
cowplot::plot_grid(
VlnPlot(test1.seur, features = check.m, group.by = "intAnno1", assay = "RNA") + NoLegend() + labs(x="Stst") +
  scale_fill_manual(values = color.A1),

VlnPlot(test2.seur, features = check.m, group.by = "intAnno1", assay = "RNA") + NoLegend() + labs(x="CR7d") +
  scale_fill_manual(values = color.A1),

VlnPlot(test1.seur, features = check.m, group.by = "intAnno1", assay = "RNA", split.by = "cnt") + NoLegend() + labs(x="Stst") +
  scale_fill_manual(values = color.test1),

VlnPlot(test2.seur, features = check.m, group.by = "intAnno1", assay = "RNA", split.by = "cnt") + NoLegend() + labs(x="CR7d") +
  scale_fill_manual(values = color.test2),
#
ncol = 2)

##
# use coord_cartesian(ylim()) instead of ylim()

cowplot::plot_grid(
VlnPlot(test1.seur, features = check.m, group.by = "cnt", cols = color.test1, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("All - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,6)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Stst.CTL","Stst.CKO")),
                               label.y = c(4.5),
                               size=3.5
                               ),

#
VlnPlot(subset(test1.seur, subset=intAnno1 %in% c("IN1")), features = check.m, group.by = "cnt", cols = color.test1, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IN1 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,6)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Stst.CTL","Stst.CKO")),
                               label.y = c(4.5),
                               size=3.5
                               ) ,

VlnPlot(subset(test1.seur, subset=intAnno2 %in% c("IN2")), features = check.m, group.by = "cnt", cols = color.test1, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IN2 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,6)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Stst.CTL","Stst.CKO")),
                               label.y = c(4.5),
                               size=3.5
                               ),

VlnPlot(subset(test1.seur, subset=intAnno2 %in% c("IMN3")), features = check.m, group.by = "cnt", cols = color.test1, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IMN3 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,6)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Stst.CTL","Stst.CKO")),
                               label.y = c(2.75),
                               size=3.5
                               ),
ncol = 2)

##
# use coord_cartesian(ylim()) instead of ylim()

cowplot::plot_grid(
VlnPlot(test2.seur, features = check.m, group.by = "cnt", cols = color.test2, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("All - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,6)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("CR7d.CTL","CR7d.CKO")),
                               label.y = c(4.5),
                               size=3.5
                               ),

#
VlnPlot(subset(test2.seur, subset=intAnno1 %in% c("IN1")), features = check.m, group.by = "cnt", cols = color.test2, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IN1 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,6)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("CR7d.CTL","CR7d.CKO")),
                               label.y = c(4.5),
                               size=3.5
                               ) ,

VlnPlot(subset(test2.seur, subset=intAnno2 %in% c("IN2")), features = check.m, group.by = "cnt", cols = color.test2, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IN2 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,6)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("CR7d.CTL","CR7d.CKO")),
                               label.y = c(4.5),
                               size=3.5
                               ),

VlnPlot(subset(test2.seur, subset=intAnno2 %in% c("IMN3")), features = check.m, group.by = "cnt", cols = color.test2, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IMN3 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,6)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("CR7d.CTL","CR7d.CKO")),
                               label.y = c(2.75),
                               size=3.5
                               ),
ncol = 2)

Nrg1

check.m <- "Nrg1"

#
cowplot::plot_grid(
VlnPlot(test1.seur, features = check.m, group.by = "intAnno1", assay = "RNA") + NoLegend() + labs(x="Stst") +
  scale_fill_manual(values = color.A1),

VlnPlot(test2.seur, features = check.m, group.by = "intAnno1", assay = "RNA") + NoLegend() + labs(x="CR7d") +
  scale_fill_manual(values = color.A1),

VlnPlot(test1.seur, features = check.m, group.by = "intAnno1", assay = "RNA", split.by = "cnt") + NoLegend() + labs(x="Stst") +
  scale_fill_manual(values = color.test1),

VlnPlot(test2.seur, features = check.m, group.by = "intAnno1", assay = "RNA", split.by = "cnt") + NoLegend() + labs(x="CR7d") +
  scale_fill_manual(values = color.test2),
#
ncol = 2)

##
# use coord_cartesian(ylim()) instead of ylim()

cowplot::plot_grid(
VlnPlot(test1.seur, features = check.m, group.by = "cnt", cols = color.test1, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("All - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,6)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Stst.CTL","Stst.CKO")),
                               label.y = c(4.5),
                               size=3.5
                               ),

#
VlnPlot(subset(test1.seur, subset=intAnno1 %in% c("IN1")), features = check.m, group.by = "cnt", cols = color.test1, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IN1 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,6)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Stst.CTL","Stst.CKO")),
                               label.y = c(4.5),
                               size=3.5
                               ) ,

VlnPlot(subset(test1.seur, subset=intAnno1 %in% c("IN2")), features = check.m, group.by = "cnt", cols = color.test1, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IN2 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,6)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Stst.CTL","Stst.CKO")),
                               label.y = c(4),
                               size=3.5
                               ),
VlnPlot(subset(test1.seur, subset=intAnno1 %in% c("IN3")), features = check.m, group.by = "cnt", cols = color.test1, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IN3 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,6)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Stst.CTL","Stst.CKO")),
                               label.y = c(4),
                               size=3.5
                               ),

VlnPlot(subset(test1.seur, subset=intAnno1 %in% c("IPAN1")), features = check.m, group.by = "cnt", cols = color.test1, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IPAN1 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,6)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Stst.CTL","Stst.CKO")),
                               label.y = c(4.5),
                               size=3.5
                               ),
VlnPlot(subset(test1.seur, subset=intAnno1 %in% c("IMN4")), features = check.m, group.by = "cnt", cols = color.test1, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IMN4 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,6)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Stst.CTL","Stst.CKO")),
                               label.y = c(3.25),
                               size=3.5
                               ),
ncol = 2)

##
# use coord_cartesian(ylim()) instead of ylim()

cowplot::plot_grid(
VlnPlot(test2.seur, features = check.m, group.by = "cnt", cols = color.test2, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("All - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,6)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("CR7d.CTL","CR7d.CKO")),
                               label.y = c(4.5),
                               size=3.5
                               ),

#
VlnPlot(subset(test2.seur, subset=intAnno1 %in% c("IN1")), features = check.m, group.by = "cnt", cols = color.test2, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IN1 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,6)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("CR7d.CTL","CR7d.CKO")),
                               label.y = c(4.5),
                               size=3.5
                               ) ,

VlnPlot(subset(test2.seur, subset=intAnno1 %in% c("IN2")), features = check.m, group.by = "cnt", cols = color.test2, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IN2 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,6)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("CR7d.CTL","CR7d.CKO")),
                               label.y = c(4),
                               size=3.5
                               ),
VlnPlot(subset(test2.seur, subset=intAnno1 %in% c("IN3")), features = check.m, group.by = "cnt", cols = color.test2, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IN3 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,6)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("CR7d.CTL","CR7d.CKO")),
                               label.y = c(4),
                               size=3.5
                               ),

VlnPlot(subset(test2.seur, subset=intAnno1 %in% c("IPAN1")), features = check.m, group.by = "cnt", cols = color.test2, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IPAN1 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,6)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("CR7d.CTL","CR7d.CKO")),
                               label.y = c(4.5),
                               size=3.5
                               ),
VlnPlot(subset(test2.seur, subset=intAnno1 %in% c("IMN4")), features = check.m, group.by = "cnt", cols = color.test2, pt.size = 0, assay = "RNA") + NoLegend() &
  geom_jitter(alpha=0.25, shape=16, width = 0.3, size = 0.25) & labs(title = paste0("IMN4 - ",check.m), x="") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,6)) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("CR7d.CTL","CR7d.CKO")),
                               label.y = c(3.25),
                               size=3.5
                               ),
ncol = 2)